On the Self‐Quenching of Relativistic Runaway Electron Avalanches Producing Terrestrial Gamma Ray Flashes

Terrestrial gamma ray flashes (TGFs) are short bursts of gamma rays occurring during thunderstorms. They are believed to be produced by relativistic runaway electron avalanches (RREAs). It is usually admitted that the number of high‐energy electrons produced in the brightest TGFs remains mostly confined within a range from 1017 to 1019. To understand the constraints in the development of RREAs, we perform self‐consistent simulations using a newly developed model with a finite acceleration region and various injection rates. We find that RREAs should naturally self‐quench for a fixed total number of runaway electrons, and hence a fixed number of bremsstrahlung photons. From the idea that TGF sources quench themselves, we derive a simple equation controlling the total number of runaway electrons. In this framework, the existence of a saturation in the electron density discovered in a previous work places a lower limit on TGF durations.


Introduction
Terrestrial gamma-ray flashes (TGFs) are short bursts of gamma-rays produced in thunderstorms by avalanches of very energetic electrons called Relativistic Runaway Electrons Avalanches (RREAs) (Dwyer et al., 2012).Several articles (e.g., Dwyer & Smith, 2005;Gjesteland et al., 2015;Mailyan et al., 2016Mailyan et al., , 2019;;Lindanger et al., 2021) showed that the amount of high-energy electrons produced during a TGF seems to consistently land between 10 17 and 10 19 .However, no explanation has been provided so far as to the physical reason for such an upper-limit.
In Gourbin and Celestin (2023b), we have reported the existence of a saturation of the low-energy electron density when the corresponding dielectric relaxation time (Maxwell time) becomes equal to the RREA characteristic time.This saturation was accompanied with a stabilization in the number of high-energy electrons and photons, suspiciously between 10 17 and 10 19 after a significant exponential growth in the RREA.In the present paper, we expand this analysis and show that the number of high-energy electrons is actually expected to be in the same range even in realistic conditions, presumably occurring in more typical TGFs.
In the next section, we present the model used and detail the method.In the third part, we present the results of the simulations.In the last part, we discuss the results and their implications on the understanding of RREAs and TGFs.

Numerical Model
The model is based on the one presented in Gourbin and Celestin (2023b).It is a self-consistent, relativistic, hybrid model using both a Monte Carlo method and an electromagnetic particle-in-cell method coupled with a fluid model solving drift-diffusion equations of low-energy charged species to simulate RREAs accurately.The particle-in-cell (PIC) method models the interaction between charged particles and the electromagnetic field.The field is evaluated on a grid from the position and motion of charged particles through a cloud-in-cell assignment scheme.The effect of the fields onto the particle follows the same interpolation scheme.The PIC method uses specific assignment and interpolation schemes to estimate the contribution of charges on the field and the effect of the field on charges.The Monte Carlo part of the code is in charge of simulating interactions between electrons and the atmosphere molecules, such as excitation of rotational, vibrational, and electronic states, elastic collisions, ionization, and attachment.For more details see (Gourbin & Celestin, 2023b) and references therein.In addition to this model, we added the treatment of high-energy photons as particles.Bremsstrahlung emission is another possible interaction for electrons, which will be the source of photons in our model.The bremsstrahlung differential cross-section is taken from Seltzer and Berger (1986).
Photons may interact with air molecules through Compton effect, photoelectric effect and electron/positron pair production.The total cross-sections used for these interactions are taken from Berger et al. (2010).The scattered photon energy and scattering angle are respectively quantified using the Klein-Nishina formula and the energy-momentum conservation law (Heitler, 1960;Lehtinen, 2000;Pilkington & Anger, 1971).For the photoelectric effect, the angle of emission is computed using Sauter's formula (e.g., Carron, 2006, Equation 2.12).As for the e + e pair production, we assume that the created electron is emitted in the same direction has the incident photon, and the positron is immediately annihilated into two 511 keV photons going in opposite directions.It is clearly a coarse approximation as positrons can travel over much longer distances in the atmosphere for the altitude range of interest here (Köhn & Ebert, 2015), but in the present work, we verified that e + e pair production occurs only on very rare occasions.

Configuration
The simulation domain in which the simulated avalanche propagates is cylindrical, axisymmetric for the electromagnetic field, with a height of 3 km and a radius of 1.2 km.The bottom part of the domain is subjected to an electric field of 16 × N N 0 kV/cm, with N 0 the density of air at the ground and N the density of air at an altitude of 12 km, the altitude we chose for all the simulations.We refer to this zone as the acceleration zone.The electric field is consistent with what is expected in streamer zones of negative leaders in positive intracloud discharges (+ICs) (e.g., Bazelyan & Raizer, 2000).The electric field is initially set to zero for z > 400 m, making the zone between z = 400 m and z = 3 km a zone where electrons interact with the air, slow down, and influence the electromagnetic field around them.Seed electrons are injected continuously at a constant rate in the acceleration zone, through a hemispherical region of radius 42 m centered on the (0,0) point (see Figure 1).Their initial energy is 1 MeV.These seed electrons can be thought of as produced by cosmic ray secondaries (e.g., Chilingarian et al., 2021), thermal runaway electrons injected by a lightning leader (e.g., Celestin et al., 2015) or relativisticfeedback-produced electrons (e.g., Dwyer, 2003;Pasko et al., 2023).The injection rate is a variable parameter between the different simulations.The different injection rates are 2.94 × 10 p , with p taking the following values: 13, 14, 15, 16, 17, 18, 19, 20, and 21.The simulations last long enough so that the avalanche reaches the end of the acceleration zone and propagates beyond that point (see Figure 1).

Results
The typical evolution of the avalanche is shown Figure 1.The injection rate used in this particular case is 2.94 × 10 21 el/s, the highest among all cases.
In all simulations, after around 1.3 μs, the avalanche reaches the end of the acceleration region.As soon as the avalanche arrives in the region where the electric field is null, it starts dispersing, as no electric field guides the electron movement anymore.Electrons still propagate over about one hundred meters before losing all their energy and eventually stop.We can also notice in this run a saturation of the low-energy (<1 keV) electron density as highlighted in Gourbin and Celestin (2023b), as well as a backward propagating field structure also discussed in the same article.This structure can be better seen in Figure 2, where the various runs stopped at different stages regarding the RREA saturation state: the cases with an injection rate below 2.94 × 10 17 electrons per second do not reach saturation before the end of the acceleration zone.The case with an injection rate of 2.94 × 10 18 el/s is close to reaching it, and starts to display the backward propagating structure, between 200 and 300 m.The three cases with higher injection rates all reach saturation before the first injected electrons leave the acceleration zone, and display different dynamics.
The number of high-energy electrons and resulting number of photons as a function of time is shown in Figure 3.In panel (a), we see the number of high-energy electrons in the simulation domain (≥1 MeV) as a function of time.
For the cases with an injection rate below 2.94 × 10 17 el/s, the number of high-energy electrons keeps increasing exponentially until the avalanche reaches the electric field cut, and the system then reaches a steady state.When they arrive at this point, the electrons at the front quickly decelerate and lose all their energy.The maximum  For the cases with the three highest injection rates, we see that the number of high-energy electrons stabilizes before the other cases, that is, before reaching the end of the acceleration zone.This stabilization is caused by the mechanism of saturation (Gourbin & Celestin, 2023b) and the associated collapse of the electric field, which constrains the number of electrons.An interesting feature is that, for these three cases, the number of high-energy electrons all lie between 10 17 and 10 18 , while the non-saturated cases all reached different magnitudes at the end of the simulations proportionally to the injection rate.
The right panel of Figure 3 shows the number of photons with an energy greater than 1 keV as a function of time for the different cases.While the curves are slightly different from the number of electrons, they follow a similar pattern: the non-saturated cases show a slowdown of the photon number at different magnitudes proportionally to the injection rate, while the saturated cases all lay between 10 17 and 10 18 photons for an increase in the injection rate by two orders of magnitude.In the following section, we argue that despite the different expected timescales and processes at play in real TGFs, a similar total number of high-energy electrons and photons to those in the saturated cases would be reached.

Discussion
In all saturated cases, the number of energetic electrons and photons stabilizes around the same order of magnitude.In our simulations, this order of magnitude is near 10 17 -10 18 .This dynamics can be better viewed on Figure 4, where the number of high-energy electrons and photons at the end of simulations are shown as a function of the number of injected electrons.
A break of the linearity of the system can be seen for cases reaching saturation (injection rates >3 × 10 17 el/s): while the number of high-energy electrons and photons at the end still slowly increases with the number of injected electrons, the increase is slowed down much more significantly, reaching for 5 × 10 15 injected electrons, 4 × 10 17 high-energy electrons at the end of the simulation.
The number of high-energy electrons can be compared with the estimates deduced from TGF observations (e.g., Mailyan et al., 2016, Figure 8;Mailyan et al., 2019, Figure 9), where the number of high-energy electrons after propagation is evaluated as mainly between 10 17 and 10 19 , with almost none beyond 10 20 .This seems to match well with the results of our saturated cases, where the electron number stagnates between 10 17 and 10 18 , especially considering that those articles use very bright TGFs (i.e., bright enough to allow spectral analysis of individual TGFs).The physical parameters used in the present study are not unrealistic.For instance, the magnitude of the electric field and the length of the acceleration region are close to those expected to be present in streamer zones of negative leaders (e.g., Bazelyan & Raizer, 2000) in positive intracloud discharges (+IC) forming an electric potential drop of ∼300 MV (Celestin et al., 2015).In such configuration, the amplification factor in the number of high-energy electrons is about 10 6 (e.g., see Figure 3), which results in the establishment of a RREA spectrum (Celestin et al., 2015).
In the present work, the magnitudes of the injection rates have been chosen to demonstrate both linear (<3 × 10 17 electrons/s) and saturation (>3 × 10 17 electrons/s) regimes.In the latter cases, RREAs lead to a quenching of the electric field over the timescale of the simulation (∼1.6 μs) while the former cases are expected to quench the electric field over longer timescales (see below).In reality the typical total duration of TGFs as observed from space is closer to ∼100 µs with rise times of a few tens of microseconds (e.g., Fishman et al., 2011;Foley et al., 2014;Lindanger et al., 2020;Xu et al., 2019).For the sake of simplicity, one can consider 15 μs as a characteristic timescale for the rise of the TGF source as derived from radio observations (e.g., Berge et al., 2022;Pu et al., 2019).
For a total number of high-energy electrons of ∼10 17 in a TGF (e.g., Dwyer & Smith, 2005), considering an amplification factor of 10 6 and a rise time of ∼15 μs, the initial injection rate can therefore be considered to be on the order of ∼10 16 electrons/s (purple case in Figures 2 and 3).It is therefore expected that in reality, TGF sources do not necessarily reach the electron saturation regime.
It is critical to consider electron attachment processes in air at an altitude of 12 km.Under an electric field of magnitude 16 N N 0 kV/cm, electron attachment is dominated by the dissociative attachment process (2-body attachment of frequency ν att2b ) with a characteristic timescale of 1 ν att ∼ 0.5 μs, even though more generally ν att = ν att2b + ν att3b , where ν att3b is the 3-body attachment frequency.Assuming an injection rate of 10 16 electrons/ s, it becomes clear that the low-energy electron density reaches a steady state at all locations in the simulation domain as an equality is reached between the ionization rate (produced by runaway electrons) and the attachment rate, this equality being reached over a duration on the order of 1 ν att ∼ 0.5 μs.Indeed, the dynamics of the lowenergy electron density n e is governed by the following equation at any given location: where F is the local z-dependent flux of high-energy runaway electrons and λ i ∼ 2 mm is the total ionization mean free path for these electrons estimated from the ionization cross-section (e.g., Celestin & Pasko, 2010;Dwyer & Babich, 2011;Kim et al., 2000).At the end of the acceleration region (z = 400 m), in the case with an injection rate of 2.94 × 10 16 s 1 , the steady state is reached for n e = F z=400m λ i ν att ∼ 10 13 m 3 (in good agreement with the results presented in Figure 2).Even though the electron density does not vary any longer, the flux of runaway electrons is still present after this time and so does the associated production of positive ions through ionization and negative ions through attachment.As a result of the flux of runaway electrons, the ion density increases linearly in time because ion recombination processes occur over a much longer timescale (see discussion in Berge et al. (2022), Section 2.1).Gourbin and Celestin (2023b) showed that the electron saturation density caused by self-consistent effects occurs when the density is high enough so that the associated relaxation time is equaling the RREA characteristic growth time: n sat e = ε 0 q e μ e ν RREA ≃ 10 15 m 3 ; ν 1 RREA = 5.16 × 10 7 s (2) where q e is the electron charge, μ e is the electron mobility, and ν RREA is the runaway electron production frequency.Figure 2 clearly shows that the quenching of the electric field is indeed rapidly observed for n e > 10 14 m 3 .In fact, the saturation density n sat e may be thought of as an electron density upper limit.In general, the electric field is quenched over the dielectric relaxation time (also named Maxwell time) driven by the local Geophysical Research Letters 10.1029/2023GL107488 conductivity, which is mostly the result of the electron density and mobility (i.e., for n ± < μ e μ ± n e , where n ± is the density of positive and negative ions and μ e μ ± ∼ 100): (3) Assuming that TGFs self-quench themselves, one can also use this timescale as a TGF characteristic rise time to estimate the corresponding maximum electron density, that is τ ∼ 15 μs.Equation 3shows that in that case, n e ∼ 10 13 m 3 , which is also consistent with the case of an injection rate of 10 16 electron/s (see Figure 2).
The number of high-energy electrons can be written as: where S is the area of the RREA at the end of the acceleration region.
As previously discussed, owing to the steady state of the electron density F z = 400m = n e ν att λ i .Substituting F and τ (Equation 3) in Equation 4, one obtains an expression for the total number of high-energy electrons in a TGF that does not depend on the low-energy density or the injection rate: For a broad range of electric fields and assuming a compact injection of runaway electrons as instantiated in the present work, this expression yields N e ∼ 10 17 .It is valid in subcritical conditions (electron density lower than saturation density) and simply comes from the assumption that the TGF quenches itself over a duration equal to the dielectric relaxation timescale caused by the production of secondary electrons.It is presumably the reason why previous works using self-consistent calculations have obtained a consistent number of TGF photons despite significant differences in the configurations, parameters, and methods employed (e.g., Dwyer, 2012;Gourbin & Celestin, 2023b;Liu & Dwyer, 2013).
The role of the conductivity increase and corresponding dielectric relaxation of the field in TGF dynamics was already explicitly mentioned and demonstrated in previous works (e.g., Dwyer, 2012).Dwyer (2007) derived a formula for the fluence, in the context of relativistic feedback, that yield similar results (see Equation 39in Dwyer ( 2007)).Equation 5 can be seen as a special case for a constant injection flux of runaway electrons, should these seed electrons be from thermal runaway, cosmic ray secondaries, or feedback processes, which is clearly instantiated by the modeling results presented in this paper.
We cannot exclude that the shortest observed TGFs (τ ≲ 10 μs) reach the electron saturation regime.In such supercritical conditions, the expected number of high-energy electrons is also appreciably constant at ∼10 17 (e.g., see Figure 3).In the present work, this would correspond to the cases with initial injection rates >2.94 ⋅ 10 18 electrons/s.This naturally indicates that TGFs cannot be shorter than ∼1 μs as determined from the RREA timescale.
The existence of a maximum electron density n sat e ∼ 10 15 m 3 reachable in a TGF and a corresponding minimum TGF timescale of ∼1 μs (RREA rate for a field of 16 × 10 5 N N 0 V/m at 12 km) are also demonstrated for the first time by the present work.
Even with more realistic TGFs consistent with maximum electron densities of n e ∼ 10 13 m 3 , we can infer a strong impact of the preionization on subsequent leaders and streamers propagation dynamics and hence the related radio observations.Because the cross-section of the RREA at the end of the avalanche is controlled by the diffusion of runaway electrons, the above analysis points to the compactness of TGF sources with radii ≲500 m.

Conclusions
Using self-consistent simulations, we highlight a mechanism limiting the number of high-energy electrons produced during a RREA and draw the following conlcusions:

10.1029/2023GL107488
• There is a maximum low-energy electron density reachable in TGF sources when RREAs reach saturation: n sat ∼ 10 15 m 3 • There exists a corresponding minimum TGF timescale equal to the RREA timescale: τ min ∼ 1 μs.
• The self-quenching of the TGF sources implies a maximum number of electrons N e ≃ ε 0 ν att λ i q e μ e ⋅ S ∼ 10 17 , weakly dependent on the electric field in the acceleration region.• In the case of typical TGFs, we can infer a strong impact of the preionization on subsequent leader/streamer dynamics following the emission of TGFs.

Figure 1 .
Figure 1.Cross-sectional view of the low-energy (<1 keV) electron density at the end of the simulation, 1.67 μs after the run started, for an injection rate of 2.94 × 10 21 s 1 .The low-energy electron density is represented by the color map.The black dots represent the location of energetic electrons at this moment of time, with an energy above 1 keV.

Figure 2 .
Figure 2. Low-energy (<1 keV) electron density profile for different injection rates at the end of the simulations (t = 1.67 μs).

Figure 3 .
Figure 3. Left: Number of high-energy (≥1 MeV) electrons as a function of time for different injection rates.Right: Number of photons (≥1 keV) as a function of time for different injection rates.

Figure 4 .
Figure 4. Number of electrons and photons at the end of the simulations as a function of the number of electrons injected in the domain throughout the duration of the simulations.