Relativistic Runaway Electron Avalanche Development Near the Electric Field Threshold in Inhomogeneous Air

Relativistic Runaway Electrons Avalanches (RREAs) development depends on the applied electric field and the environment's air density. This dependency controls the RREA exponential growth length scale. The RREA development affects the bremsstrahlung excess occurring due to the passage of charged particles through the thundercloud's electric fields, the gamma‐ray glow. We used Monte Carlo simulations to develop an empirical model showing the RREA behavior in a realistic atmospheric density profile. The new formulation shows how the density variation modulates the electron population under electric field strengths near the RREA electric field threshold. The model limits the initial RREA altitude range as a function of the electric field strength. The new model is valid between ∼0.6 and ∼18 km, covering the relevant heights to investigate the generation of ground‐detected gamma‐ray glows.

Such E S makes the RREA development sensitive to air density variations because the avalanche have an initial altitude, h i , where E S > E th (h i ) and travels to regions where E S < E th (h i ).Thus, inhomogeneous air will change both the RREA and the feedback mechanism.We present E S normalized by E th0 in our analysis for convenience, that is, δ E (E S ) = E S /E th0 .Figure 1 shows how the density profile creates vertical regions where the RREA is allowed due to altitude dependency with the vertical distance from source Δh = h i − z.
The current work investigates RREA development in a realistic atmospheric vertical density profile.It shows the transition between RREA growth and decay associated with the electric field strength, E S , and the altitude z, see Figure 1.Diniz et al. (2022) showed the extension of electrons spatial range in ambients where E S < E th (z), in which there is no avalanche but still a longer-lasting high-energy electron population.Now, we evaluate the regime E S > E th (h i ) transitioning to E S < E th (z) due to the vertical density change.This transition is relevant to gamma-ray glow observations during Japanese Winter (Wada et al., 2021) modeled by Diniz et al. (2023).We show in Diniz et al. ( 2023) that E S ∼ E th (z) are suitable for the Gamma-Ray Observations of Winter THunderclouds (GROWTH) reports (Wada et al., 2021).Here, we present the third of our paper series pursuing a theoretical framework for gamma-ray glows.The first one is Diniz et al. (2022), and the second paper is Diniz et al. (2023).They are referred here to as paper I and II, respectively.We describe the number of relativistic electrons,  e− , as a function of E S , h i , and Δh; thus,  e− = e− (Δℎ, ℎi, S) .Finally, we show the physical consequences of the inhomogeneous air density upon RREA and discuss the corollaries regarding gamma-ray glow generation.

Simulation Setup and Analysis Method
Our simulations use GEometry ANd Tracking 4 (GEANT4) version 10.4.3 with the standard "FTFP_BERT_ EMZ" physics list (Agostinelli et al., 2003;Allison et al., 2006Allison et al., , 2016) ) including all the relevant physics.The simulation geometry is a cylinder with a 15 km radius and different heights determined by the number of air layers with different discrete densities following the Naval Research Laboratory Mass Spectrometer Incoherent Scatter radar (NRLMSIS) atmospheric model (https://kauai.ccmc.gsfc.nasa.gov/instantrun/nrlmsis/). We used the standard air composition, which consists of 78.085% nitrogen, 20.950% oxygen, and 0.965% argon.The geometry includes a vertical, upwards homogeneous electric field filling the ambient varying the strength E S between 0.220 and 0.270 MV/m in 0.005 MV/m steps.The electric field is homogeneous and time-constant because existing field measurements in thunderclouds do not provide sufficient data of an accurate  ⃗  profile for our simulation.The air layers have a standard 50 m thickness of different densities.The initial altitude of the point pencil-like electron beam, h i , varies from 0.75 to 2.50 km in 0.25 km steps.The initial particles are fifty 20 MeV electrons to ensure the avalanche development.The electrons recording are at each interface between two air layers.The simulations include a 100 keV internal energy cutoff for all particles except positrons due to annihilation possibility.We adjust the internal energy cutoff according to our data grid to reduce computational costs.The maximum altitude considered in our data is 2.5 km with an air density of ∼0.95 kg/m 3 .Electrons with 100 keV kinetic energy cannot seed RREA at this density under a 0.270 MV/m electric field applied, especially with the increasing density ambient.We count the electrons with energy above 1 MeV, which can contribute to the RREA in the electric field, E S ∼ E th (z) (Sarria et al., 2018).Finally, we normalize the count by the initial electron number.
Generally, the relativistic electron number,  e− , depends on the position,    , the electric field  ⃗  and, the electron kinetic energy ɛ.The variation with energy can be neglected in the RREA population dynamics since the relativistic electron multiplication process is the Møller scattering.Its cross-section is inversely proportional to the primary particle velocity, , where c is the speed of light and v 1 the primary electron velocity making the production of relativistic electrons constant as the incoming electron energy increases (Sarria et al., 2018).Finally, we approximate the position dependency to the vertical direction, z, aligned with the applied electric field.

Electron Avalanche Development as a Function of E S and h i
The RREA population growth in a homogeneous medium is exponential,  e− ∝  ∕ R (S,r) , where λ R (E S , n r ) is the characteristic avalanche length defined in Equation 1 in meters (Dwyer, 2003) with (1) The atmospheric density profile has a scale height of approximately 8 km (Picone et al., 2002).Thus, altitude variations of 1 km correspond to ∼10% change in density.The homogeneous density approximation is reasonable if E S ≫ 0.276n r according to Equation 1.However, the density variations are relevant if E S ∼ E th (z) because the electron population transitions between RREA region and damping regions, see Figure 1.Sarria et al. (2023), and paper II relate the latter electric field regime to gamma-ray glows.Thus, understanding the RREA dynamics in an inhomogeneous air density profile is necessary to describe the gamma-ray glow source.
Figure 2 shows the simulated number of electrons with kinetic energy higher than 1 MeV,  e− (Δℎ, ℎi, S) , versus the vertical distance from the initial altitude, Δh, the legend indicates the E S and h i in each case.Panels A refer to  e− (Δℎ, 2 km, S) while panels B display  e− (Δℎ, ℎi, 0.265 MV∕m) .
Equation 2 fits Figure 2 data.The fitting becomes more accurate for increasing E S due to a better avalanche development for E S increasingly higher than (2) In Equation 2, λ G (E S , h i ) and λ D (E S , h i ) are the growth length and damping length, respectively.The fitting constants determine the values of λ G (E S , h i ) and λ D (E S , h i ) for each simulation.
In particular, the RREA threshold at 2 km is ∼0.231MV/m while the avalanche development starts to be noticeable for E S ∼ 0.240 MV/m, in our data set.Stronger electric fields allow RREA development in shorter space.
Likewise, the fit quality also improves for increasing h i , for the same reason.An increasing difference between E S and E th (h i ) provides a wider space for avalanche development.The electron shower enters the damping region before an avalanche development once E S ≈ E th (h i ) because of increasing density, see Figure 1.
We obtain the altitude of maximum  e− (Δℎ, ℎi, S) , z max (E S , h i ) by differentiating Equation 2 with respect to Δh, resulting in Equation 3, (3) Thus, the electron number,  e− , variation with the altitude Δh can be written as, and the adapted avalanche characteristic length, Λ(Δh, h i , E S ), is, Figure 3 shows the dependence of λ G (E S , h i ), λ D (E S , h i ), and Δh max (E S , h i ) behavior as a function E S and h i , separately.Following Equation 2, λ G (E S , h i ) must approach λ R (E S , n r ) for either large E S and large h i to recover the exponential growth previously established (Dwyer, 2003) because these conditions induces the regime E S ≫ E th (h i ).Combining the electric field strength analysis of Figure 3A with the initial shower altitude evaluation of Figure 3B we reach the multi-variable forms, In Equations 6 and 7, A i and B i , with the subindex "i" relative to either growth (G) or damping (D) quantities, are the multi-variable fitting constants.A two variable fitting within our valid data grid (E S > E th (h i ) regime) provides, A G = 0.014 ± 0.013, B G = 0.256 ± 0.030 km, A D = 1.508 ± 0.241, B D = 0.760 ± 0.422 km.

Discussion and Conclusion
The atmospheric density variation damps the RREA development for electric field strengths close to E th (h i ).The RREA region is smaller than the electric field region because the electron multiplication starts at an altitude of E S > E th (h i ) and descend to locations where E S < E th (z).An analog dynamic should be present throughout the thundercloud complex electrical structures such as the thundercloud reactor model (Stadnichuk et al., 2021).
Equations 6 and 7 imposes geometrical limitations upon Equation 2 model as λ G (E S , h i ) and λ D (E S , h i ) need to be positive to have physical meaning.These geometrical limitations are, H S (E S ) and H I (E S ) are the superior and inferior limits for the initial shower altitude, H S (E S ) ≥ h i ≥ H I (E S ).The limits H S (E S ) and H I (E S ) further constrain the RREA region.Moreover, the range of possible RREA initial altitudes reduces as weaker is the electric field because  lim S →0S(S) = 0 and  lim S →0I(S) = ∞ .Considering the parameter values found at the end of Section 3, the boundaries for the initial shower altitude from Equations 8 and 9 become H S (E S ) = (18.286± 17.114)δ E (E S ) km and H I (E S ) = (0.504 ± 0.291)/δ E (E S ) km.The errors result from uncertainty propagation from the fitting parameters.These values are relative to our discrete data grid of E S and h i .Thus, a denser data grid will reduce the errors and change the boundaries.In the context of this work, our main interest is in gamma-ray glow ground measurements correlated with winter Japanese thunderclouds (Wada et al., 2021).Hence, the boundary values cover our purposes.
Figure 4 evaluates the model boundaries and validity.Panel A shows z max (E S , h i ) from Equation 3 as a function of E S and h i .The conditions that limit the model applicability, E th (h i ), H S (E S ), and H I (E S ) are overplotted as dashed using Equation 3for the data points that fulfill the condition E S > E th (h i ).
There are 88 simulated configurations, 52 satisfying the validity condition.Of the 52 valid ones, 69.2% match the fitting by a 10% difference and 80.7% within a 20% difference.Panel B.1 outliers points show that the data and fitting agreement is sensitive to how stronger E S is than E th (h i ).Panel B.2 shows the normalized difference between the used E S and E th (h i ).The electric field strength, E S , must be higher than E th (h i ) to assure the model applicability, thus, all the configurations above zero are valid.
There are two error sources: (a) the model does not contemplate the particles' horizontal displacement, and (b) the two-dimensional fitting leads to relatively large error bars in comparison to the average values due to the valid data grid shape, see Figure 4 panel (B.2).The error (1) source is negligible as the RREA threshold (E th (z)) is a result of the horizontal displacement Dwyer (2003).Error (2) makes the single dimensional fittings of λ G (E S , h i ) and λ D (E S , h i ) (Figure 3) more accurate than the two dimensional fitting (Figure 4).Nevertheless, we show the two-dimensional fitting to visualize the model corollary, z max (E S , h i ).
The model validity up to 18 km embraces the relevant altitudes regarding thundercloud heights.It evaluates the RREA behavior with non-uniform density ambient.The green curve, H S (E S ), of Figure 4 reach unrealistic values for the uniform electric field extension hypothesis since it covers lengths beyond any acceleration region.Thus, all realistic initial altitude, h i , are below H S (E S ).Nevertheless, we show the full altitude range to display the model's validity completely.
Equations 3 and 7 show the z max (E S , h i ) behavior regarding the electric field strength and the initial altitude, max (𝐸𝐸S, ℎi) ℎi Equation 10 shows a global maximum as a function of the electric field strength, E S , where E S = 2E th0 /A D .
Considering A D = 1.508 ± 0.241, Equation 10 is positive for E S below ∼ 0.367 ± 0.060 MV/m.Thus, z max (E S , h i ) will increase until E S is enough to well develop the RREA.The fully developed RREA will decrease z max (E S , h i ) because the electrons are moving downwards and the population maximum will always be the avalanche final point.On the other hand, Equation 11 is always negative regardless of h i value.Thus, an horizontal line in Figure 4 panel (A) will present a negative z max (E S , h i ) rate that decreases with larger h i .
There is some freedom on where to consider the initial shower altitude since there is a constant cosmic-ray flux hitting thunderclouds at all times.The initial shower altitude is at the beginning of the electrified space.In the gamma-ray glow framework, This point would be the top of the region between the negative charge center and the lower positive charge center at the cloud base considering the classical tripolar cloud structure (Takahashi, 1978;Williams, 1989) for ground detection such as (Wada et al., 2021).
Winter thunderclouds in Japan have base altitude of 0.2-0.8km (Takahashi, 1978) which allowed several gammaray glow observations (Torii et al., 2002;Tsuchiya et al., 2012;Wada et al., 2019Wada et al., , 2021)).Recently, Paper II reported geometric requirements to explain the gamma-ray glow spectra, reported by Wada et al. (2021), considering air with a homogeneous density of ground level.The configurations include three geometric coordinates, acceleration region of vertical length, H E , electric field strength, E S and, an attenuation region of vertical length, H a , making a triad of (E S , H E , H a ).The coordinates triad that produces spectra closest to the observations are (0.31 MV/m, 1 km, 0.4 km) and (0.30 MV/m, 0.9 km, 0.3 km).We must transpose these coordinates to the inhomogeneous air situation.The real electric field strength, E S , may be weaker than proposed in Paper II due to lower air density in source altitudes.The electric field length, H E , should also be modified to embrace realistic potentials.In particular, the conjugate variables (E S , H E ) will determine together the feedback mechanism possibility.
The feedback mechanism depends on the avalanche length as in Equation 1 while the adapted avalanche length is Equation 5 in inhomogeneous air.This modification alters the feedback dynamics.The dependency of Equation 5on the altitude makes the exponential decreases for z < z max .The maximum  e− (Δℎ, ℎi, S) altitude, max , marks the transition to between RREA and damping regions.The electric field is less capable of accelerating electrons downward and accelerating positrons upwards as the RREA descends to higher-density altitudes, that is, the damping region.This transition will then decrease the probability of a seed electron producing second-generation seeds near the start of the avalanche region as described in (Dwyer, 2007).The increasing density increases the friction, leaving less energy available for backscattered X-rays.Dwyer (2007) states that the general principles of the feedback mechanism are the same regardless of the electric field configuration.This principle also holds for the density variation while the probabilities of backscattering for the feedback change substantially.In realistic cases with smaller electric field regions, the density variation holds the feedback mechanism down.It should rely on electric fields associated with intense gamma-ray glows.Further development on the feedback mechanism in this context requires a follow-up dedicated analysis of future works exploring the electric field size and strength.
Both evaluations, paper II and the present work, consider homogeneous electric fields within the acceleration region.The electric field may vary in time and space.Such variations would modify the RREA dynamic in reality.Thus, the missing points to resolve the gamma-ray glow process are the electric field accurate geometry and the electric field variation with time inside the thundercloud.

Figure 1 .
Figure 1.E th (z) represented as a function of the vertical distance from the source, Δh, (black line) compared with different electric field strength, E S , values represented by the colored dashed lines and translated into electric field strength normalized by the Relativistic Runaway Electrons Avalanche threshold at ground level (E th0 ), δ E (E S ).The gray shaded area marks the damping region.This figure covers the initial altitude, h i , range of 0-2.5 km.

Figure 3
Figure 3 panels (A.3) and (B.3) display the model accuracy through comparison between the Δh max from fitting  ( Δℎ fit max ) and data  ( Δℎ data max ), both normalized by h i .The convergence of Δh max to h i shows the RREA better development with an increasing difference between E S and E th (h i ).Resulting in the approximation of Equation 2 to the already establishedDwyer [2003] empirical formulation for high E S or h i , see Figure3panels (A.1) and (B.1).

Figure 2 .
Figure 2. Number of electrons with kinetic energy higher than 1 MeV,  e− (Δℎ, ℎi, S) , as a function of vertical distance from source, Δh, for different electric field strength, E S , in panel (A.1) and different initial altitude, h i , in panel (B.1)-each one is indicated in the legends.The dots are simulation results, and the curves are fitting models by Equation 2. The error bars follow Poisson statistics as in papers I and II.Panels (A.2) and (B.2) show the respective data-to-fit ratio.

Figure 3 .
Figure 3. Growth (λ G (E S , h i )) and damping (λ D (E S , h i )) lengths, and the vertical distance between the  e− (Δℎ, ℎi, S) maximum and the initial altitude, h i , Δh max (normalized by h i ), as functions of E S and h i for the cases in Figure 2. Panels (A.1) and (B.1) compare the respective growth lengths with the characteristic avalanche length, λ R (E S , n r ).Panels (A.2) and (B.2) show the damping (λ D (E S , h i )) lengths.Panels (A.3) and (B.3) compare Δh max retrieved from the data  ( Δℎ data max

Figure 4 .
Figure 4. Panel (A) shows the height of avalanche maxima z max (E S , h i ) as a function of electric field strength E S and altitude z, displayed by the contour lines with values in kilometers.The dashed red line marks the Relativistic Runaway Electrons Avalanche avalanche threshold, E th (z).The green and blue dashed lines are the model applicability superior and inferior limits H S (E S ) and H I (E S ) showing the valid region for the present model, respectively.Panel (B) displays the model validity regarding our discrete data grid.The upper panel (B.1) shows the relative difference between Δh max from the fitting model  ( Δℎ fit max ) and from the simulated data  ( Δℎ data max