Relativistic runaway electron avalanches: unified density-dependent scaling and transport

Relativistic runaway electron avalanches (RREAs) are commonly characterized by an empirical relation linking the avalanche length to the electric-field excess above the runaway threshold. In this study, CORSIKA simulations are used to examine the behavior of this parameterization under different atmospheric densities and electric-field conditions and to investigate the subsequent propagation of particles beyond the accelerating field region. Simulations are performed for four high-altitude observational sites.

The effective avalanche lengths derived from the vertical electron profiles exhibit inter-station scatter when fitted with a single coefficient across the combined dataset. Incorporating an additional density-dependent term reduces this scatter and improves the fit, increasing the coefficient of determination from R²≈0.90 to R²≈0.99. Simulations of particle propagation beyond the electric-field region show different attenuation behavior for the electron and gamma-ray components. Reassessment of the empirical free path distance (FPD) formulation for a 100 m field-to-detector separation yields an electron–gamma energy coefficient of C₁≈1.37, while the density dependence of electron energy losses is accounted for using stopping-power data from the National Institute of Standards and Technology (NIST) ESTAR (Stopping Powers and Range Tables for Electrons) database.

1 Introduction

High-energy particle phenomena in the atmosphere constitute a rapidly developing field at the intersection of cosmic-ray physics, atmospheric electricity, and meteorology. Over the past decades, a wide variety of high-energy atmospheric phenomena have been observed, including thunderstorm ground enhancements (TGEs), gamma-ray glows, terrestrial gamma-ray flashes (TGFs), and the recently identified flickering gamma-ray flashes (FGFs) (Chilingarian et al., 2011; Fishman et al., 1994; Dwyer et al., 2012; Østgaard et al., 2024; Marisaldi et al., 2024). Despite the diversity of detection platforms and energy ranges, mounting experimental evidence suggests that these phenomena are closely related and arise from electron acceleration in atmospheric electric fields (Dwyer and Uman, 2014). Although the microphysics of gamma-ray production, electron acceleration, and electromagnetic interactions in air are well established (Bethe and Heitler, 1934; Koch and Motz, 1959), the atmospheric environment in which these processes occur is highly dynamic and difficult to characterize, particularly because reliable measurements of the spatial and temporal structure of electric fields within thunderclouds remain limited. RREA development has been investigated using analytical, numerical, and Monte Carlo approaches (Babich et al., 2001; Dwyer, 2003; Coleman and Dwyer, 2006; Dwyer et al., 2012; Skeltved et al., 2014), while empirical and semi-empirical formulations remain useful for describing the main characteristics of energetic particle populations under atmospheric electric fields.

A key step in the theoretical description of particle acceleration under thunderstorm conditions was the development of the relativistic runaway electron avalanche (RREA) theory (Gurevich et al., 1992; Dwyer, 2003). In sufficiently strong atmospheric electric fields (AEFs), energetic electrons can overcome ionization losses and enter a runaway regime, where they continue to gain energy from the electric field and produce secondary electrons through collisions with air molecules. This process leads to an avalanche multiplication of electrons and the accompanying production of high-energy bremsstrahlung gamma rays. In addition to this conventional avalanche development, more complex RREA models include the relativistic feedback mechanism. In this mechanism, backward-propagating positrons and high-energy photons generated within an avalanche can produce new seed electrons and initiate successive avalanches. Under sufficiently strong and extended electric-field conditions, this feedback can become self-sustaining and limit the maximum electric field that can persist in air (Dwyer, 2003, 2007). A central parameter characterizing this avalanche development is the e-folding length λ (E, ρ), which describes the characteristic distance over which the number of runaway electrons grows exponentially. For the conventional RREA development considered in this study, a widely used empirical description of this process was introduced by Dwyer (2003), who proposed the relation

$$\begin{array}{}\text{(1)}& \mathit{\lambda }(E,\mathit{\rho })={\displaystyle \frac{K}{E_z-E_{\mathrm{th}}\left(h\right)}},\end{array}$$

where E_z and $E_{\mathrm{th}}\left(h\right)=E_{\mathrm{0}}\cdot \mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}$ (the runaway threshold electric field) are expressed in kV m⁻¹, λ is in meters, and K in kV. Here, ρ₀=1.225 kg m⁻³ is the air density at sea level, and ρ is the air density at altitude z. Dwyer (2003) reported $K\approx \mathrm{7.2}\times {\mathrm{10}}^{\mathrm{3}}$ kV, while subsequent Monte Carlo calculations by Coleman and Dwyer (2006) obtained a closely related parameterization for electric fields above 300 kV m⁻¹ and showed that a different empirical fit provides a better description close to the runaway threshold. Coleman and Dwyer (2006) also explicitly accounted for atmospheric density by scaling the electric field and avalanche length to their equivalent sea-level values. These results indicate that the applicability of a given avalanche-length parameterization should be considered in relation to the electric-field range and atmospheric conditions for which it is evaluated.

The concept of free path distance (FPD) has been introduced as a phenomenological quantity that characterizes the ability of energetic particles to propagate beyond the electric-field region and reach ground-based detectors (Chilingarian et al., 2021, 2022). Although this approach provides a useful link between simulations and observations, its empirical coefficients are constrained by limited datasets and simplified assumptions. In this work, two empirical relations that describe complementary stages of RREA evolution are revisited and optimized: avalanche development within the electric field and particle propagation beyond it. Both are evaluated and optimized using CORSIKA simulations for multiple high-altitude stations with different atmospheric densities. In particular, a proportionality between electron and gamma-ray energies at the field boundary is identified, and density-dependent energy losses during propagation are explicitly accounted for.

2 Simulation framework

To investigate the development of relativistic runaway electron avalanches (RREAs) and the subsequent propagation of energetic particles, Monte Carlo simulations were performed using CORSIKA (COsmic Ray SImulations for KAscade) code (Heck et al., 1998), version 7.7500, which includes the effect of atmospheric electric fields on particle transport (Buitink et al., 2009).

As in Chilingarian et al. (2026), the atmospheric electric field was modeled as a vertically oriented uniform layer with a total thickness of 2000 m above the detector level. As shown in Sect. 3, a layer thickness of 2000 m is greater than the effective avalanche lengths obtained under all simulated electric-field conditions, allowing the avalanche multiplication process to develop before particles leave the acceleration region. To resolve the vertical development of electron and gamma-ray populations, the field region was subdivided into 20 layers of 100 m each, and particle numbers were recorded at each boundary. In addition to the observation levels within the electric-field region, an additional level was defined at the lower boundary of the electric field (exit), as well as a detector level (det) located 100 m below this boundary, corresponding to the station altitude above sea level. These levels are used for the post-field transport and FPD analyses.

Simulations were performed for four high-altitude observational stations with distinct atmospheric densities, all of which have reported thunderstorm-related particle enhancements: Aragats (ρ=0.77 kg m⁻³, roughly 3000 m a.s.l. – above sea level) (Chilingarian et al., 2022, 2024); Nor Amberd (ρ=0.89 kg m⁻³, about 2000 m a.s.l.) (Chilingarian et al., 2025); Lomnický Štít (ρ=0.83 kg m⁻³, approximately 2600 m a.s.l.) (Chum et al., 2020; Kisvárdai et al., 2025); and LHAASO (ρ=0.66 kg m⁻³, around 4400 m a.s.l.) (Aharonian et al., 2023). Across the four configurations, the electric-field regions span altitudes of approximately 4000–6400 m.

The electric-field ranges used in the simulations were selected separately for each observation site to obtain sustained supercritical RREA multiplication within the 2000 m field region under the corresponding atmospheric conditions. This selection was guided by the previous CORSIKA study (Chilingarian et al., 2026), in which the dependence of the RREA initiation threshold on altitude and atmospheric density was investigated. The local runaway threshold is defined as $E_{\mathrm{th}}\left(h\right)=E_{\mathrm{0}}\cdot n\left(h\right)$, where $n\left(h\right)=\mathit{\rho }\left(h\right)/{\mathit{\rho }}_{\mathrm{0}}$ is the normalized atmospheric density and E₀=284 kV m⁻¹ is the adopted sea-level threshold, consistent with previous studies. For the empirical scaling analysis, a single representative atmospheric density is assigned to each station, corresponding to the midpoint altitude of the 2000 m electric-field layer. This reference density is used to calculate the station-specific values of n and the corresponding runaway threshold E_th for each station. Because E_th decreases with decreasing air density, the same absolute electric field does not correspond to the same degree of supercriticality at different altitudes. Therefore, the absolute applied-field values were not kept identical across the four sites but were varied taking into account the local atmospheric density and the corresponding runaway threshold. Accordingly, field ranges of 180–210 kV m⁻¹ for LHAASO, 200–230 kV m⁻¹ for Aragats, 210–240 kV m⁻¹ for Lomnický Štít, and 230–260 kV m⁻¹ for Nor Amberd were used, with the electric field varied in steps of 10 kV m⁻¹ within each range. To express the selected fields under different atmospheric-density conditions on a common physical basis, the applied field is normalized to the local runaway threshold, $E/E_{\mathrm{th}}$. This dimensionless ratio characterizes the degree of field supercriticality under the local atmospheric conditions and provides a common basis for comparing the field regimes selected for the different observation sites.

Each simulation corresponds to the injection of a single seed electron into the electric-field region. The initial energy spectrum was adopted from the EXPACS model (Sato, 2015) and follows a power-law distribution with an index of 1.173 in the energy range 1–300 MeV. To ensure statistical stability, 500–1000 simulation events were performed for each field configuration. During particle propagation, electrons and gamma rays were tracked until their energies decreased to 0.05 MeV. At each observation level within the electric-field region, the numbers of electrons and gamma rays were recorded. The electron population N_e(d) is used to determine the effective avalanche length, while particle spectra at the field boundary and detector level are used for the FPD analysis.

3 Avalanche development within the AEF

The vertical profiles of the electron population from CORSIKA simulations provide a basis for determining the effective avalanche length of relativistic runaway electron avalanches. Particle numbers were recorded at 100 m intervals within the L=2000 m electric-field layer. To minimize sensitivity to the field boundaries, the effective avalanche-length analysis was restricted to the common fitting interval d=200–2000 m. For each value of the applied electric field and each observational site, the simulations produce a depth-dependent electron population, N_e(d), where d denotes the downward propagation distance from the top of the electric-field layer. The absolute height H is then given by:

$$\begin{array}{}\text{(2)}& H=H_{\mathrm{st}}+\mathrm{2000}-d,\end{array}$$

where H_st is the altitude of the observational station above sea level. When the electric field exceeds the runaway threshold field E_th(h), the number of runaway electrons increases approximately exponentially with propagation distance. In this regime, the development of the avalanche can be described by the relation

$$\begin{array}{}\text{(3)}& N_{\mathrm{e}}\left(d\right)=N_{\mathrm{0}}{e}^{d/\mathit{\lambda }},\end{array}$$

where N₀ is the electron population extrapolated to d=0, and λ is the characteristic avalanche length. Taking the natural logarithm of the expression above yields a linear dependence of the form

$$\begin{array}{}\text{(4)}& \ln N_{\mathrm{e}}\left(d\right)=c+bd,\end{array}$$

where c=ln N₀ and $b=\mathrm{1}/\mathit{\lambda }$. Thus, the avalanche length can be obtained from the slope of the linear fit to the logarithmic electron profile.

For each simulated electric-field strength, the electron population N_e(d) was averaged over the ensemble of simulation events. The natural logarithm of the averaged electron counts was then fitted as a function of depth using linear regression. The effective avalanche length was then obtained as

$$\begin{array}{}\text{(5)}& {\mathit{\lambda }}_{\mathrm{eff}}=\mathrm{1}/b,\end{array}$$

The standard error of the fitted slope b was propagated through Eq. (5) to obtain the corresponding fit-derived uncertainty σ_λ in λ_eff.

Repeating this process for various electric-field strengths and across all observational sites produces a set of effective avalanche lengths, λ_eff, corresponding to different values of the electric-field excess, E_z−E_th(h). These values serve as the basis for optimizing the empirical coefficient K in the avalanche-length relation (Eq. 1). The extracted avalanche lengths are then compared with the dependence predicted by Eq. (1).

The CORSIKA simulations are used to examine how consistently a single coefficient K describes the effective avalanche lengths obtained for the four observational sites.

For each simulated electric-field value and for each observational site, the variable

$$\begin{array}{}\text{(6)}& x=\mathrm{1}/\left(E_z-E_{\mathrm{th}}\left(h\right)\right),\end{array}$$

was computed․ According to Eq. (1), the avalanche length depends linearly on this variable:

$$\begin{array}{}\text{(7)}& \mathit{\lambda }=Kx\mathrm{․}\end{array}$$

Thus, the coefficient K can be determined by fitting a straight line through the origin in the (x,λ_eff) parameter space.

The fitting procedure was performed using the least-squares method constrained to pass through the origin. In this case, the optimal value of the coefficient is given by

$$\begin{array}{}\text{(8)}& K=\sum _i{\mathit{\lambda }}_i{x}_i/\sum _i{x}_i^{\mathrm{2}},\end{array}$$

where $x_i=\mathrm{1}/(Ez{,}_i-E_{\mathrm{th}}(h\left)\right)$, λ_i is the corresponding effective avalanche length, and the index i labels the four simulated supercritical electric-field values used for each observational station in the least-squares fitting.

The fit-derived uncertainties of the individual effective avalanche lengths were propagated to the station-specific coefficient K. Treating the x_i values as fixed and the uncertainties of the fitted λ_i values as independent, the uncertainty in K was calculated as

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_K={\displaystyle \frac{\sqrt{\sum _i{x}_i^{\mathrm{2}}{\mathit{\sigma }}_{{\mathit{\lambda }}_i}^{\mathrm{2}}}}{\sum _i{x}_i^{\mathrm{2}}}}.\end{array}$$

Here, ${\mathit{\sigma }}_{{\mathit{\lambda }}_i}$ is the fit-derived uncertainty of the effective avalanche length λ_i for the ith electric-field configuration, x_i is defined by Eq. (6), and σ_K is the resulting propagated uncertainty of the station-specific coefficient K.

Applying this procedure consistently to the CORSIKA simulation results for all four observational sites yields the following station-specific coefficients, including the propagated fit-derived uncertainties:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{K}_{\mathrm{Aragats}}=(\mathrm{10.2}\pm \mathrm{0.4})\times {\mathrm{10}}^{\mathrm{3}}\mathrm{kV}\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{LHAASO}}=(\mathrm{11.2}\pm \mathrm{0.3})\times {\mathrm{10}}^{\mathrm{3}}\mathrm{kV}.\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{Lomnick}\mathrm{ý}\mathrm{Š}\mathrm{t}\mathrm{í}\mathrm{t}}=(\mathrm{8.9}\pm \mathrm{0.4})\times {\mathrm{10}}^{\mathrm{3}}\mathrm{kV}\mathrm{․}\\ {\displaystyle }& {\displaystyle }{K}_{\mathrm{NorAmberd}}=(\mathrm{8.6}\pm \mathrm{0.2})\times {\mathrm{10}}^{\mathrm{3}}\mathrm{kV}.\end{array}$$

The station-specific coefficients obtained from the CORSIKA simulations are larger than the $K=\mathrm{7.2}\times {\mathrm{10}}^{\mathrm{3}}$ kV value reported by Dwyer (2003), ranging from approximately 8.6×10³ to 11.2×10³ kV across the four sites. Accounting for the propagated fit-derived uncertainties does not remove the station-to-station variation; the corresponding uncertainty intervals do not collapse onto a common range, although partial overlap occurs for some neighboring values. A simultaneous fit of Eq. (1) to all 16 simulation points yields a best-fit coefficient of $K_{\mathrm{fit}}\approx \mathrm{9.58}\times {\mathrm{10}}^{\mathrm{3}}$ kV. This global fit is used below as the reference for comparison with the density-dependent formulation. For comparison with previous RREA calculations, the present avalanche lengths can also be expressed in the conventional density-scaled form, using $E_z/n$ and n⋅λ_eff. Earlier numerical studies have reported broadly comparable avalanche-length scales in these similarity-scaled variables, although quantitative differences among individual models remain, particularly near the runaway threshold (Coleman and Dwyer, 2006; Dwyer et al., 2012; Skeltved et al., 2014). As an illustrative example, for the Aragats station at E_z=210 kV m⁻¹, the present simulation corresponds to $E_z/n\approx \mathrm{334}$ kV m⁻¹ and $n\cdot {\mathit{\lambda }}_{\mathrm{eff}}\approx \mathrm{201}$ m, which is within the range of previous Monte Carlo estimates at comparable density-scaled electric fields.

The observational sites considered in this work span a substantial range of atmospheric densities, from approximately 0.66 kg m⁻³ at LHAASO to 0.89 kg m⁻³ at Nor Amberd. Although the local runaway threshold E_th(h) accounts for density through the threshold scaling, the station-to-station variation of the fitted K values persists after the fit-derived uncertainties are taken into account. This motivates testing an additional explicit density dependence in the empirical parameterization. An additional density dependence is introduced through the generalized empirical relation

$$\begin{array}{}\text{(10)}& \mathit{\lambda }(E,h)=Kn(h{)}^{\mathrm{a}}/\left(E_z-E_{\mathrm{th}}\left(h\right)\right),\end{array}$$

where $n\left(h\right)=\mathit{\rho }\left(h\right)/{\mathit{\rho }}_{\mathrm{0}}$ is the normalized air density, ρ(h) is the local atmospheric density at altitude h, ρ₀=1.225 kg m⁻³ is the air density at sea level, and a is an empirical exponent describing the additional density dependence.

Equation (10) was fitted jointly to the effective avalanche lengths from all four stations, with K and a as free parameters. The joint fit yields $K\approx (\mathrm{6.39}\pm \mathrm{0.26})\times {\mathrm{10}}^{\mathrm{3}}$ kV and $a\approx -\mathrm{0.923}\pm \mathrm{0.083}$, where the quoted uncertainties are derived from the joint fit․

The fitted exponent $a=-\mathrm{0.923}\pm \mathrm{0.083}$ indicates an additional approximately inverse dependence on atmospheric density beyond that already included through the density-scaled runaway threshold $E_{\mathrm{th}}=E_{\mathrm{0}}\cdot n$. Thus, the factor n^a represents an additional density-dependent correction within the present empirical parameterization. When ${\mathit{\lambda }}_{\mathrm{eff}}/n(h{)}^{\mathrm{a}}$ is plotted as a function of the variable $\mathrm{1}/(E_z-E_{\mathrm{th}}(h\left)\right)$, the results from all four stations show reduced inter-station scatter around an approximately linear relation (Fig. 1). Figure 1 summarizes the simulation results obtained for the four observational sites considered in this work and compares the scaling relations given by Eqs. (1) and (10).

Figure 1 Comparison of the standard and density-dependent parameterizations of the RREA avalanche length. (a) Standard parameterization for the combined four-station dataset, showing the global fit and the reference relation. (b) Density-dependent parameterization for the same dataset, showing reduced inter-station scatter after accounting for the additional density dependence.

Figure 1a shows the effective avalanche length λ_eff as a function of $\mathrm{1}/(E_z-E_{\mathrm{th}}(h\left)\right)$. The results exhibit an approximately linear trend with inter-station scatter. Figure 1b shows the same data after applying the density-dependent normalization n(h)^a. With the fitted parameters $K\approx \mathrm{6.39}\times {\mathrm{10}}^{\mathrm{3}}$ kV and $a=-\mathrm{0.923}$, the inter-station scatter is substantially reduced.

The station-specific coefficients K_station and the effective coefficients K_eff(h) are summarized in Table 1.

Figure 2 Parity comparison between the effective avalanche lengths obtained from the CORSIKA simulations and those predicted by the global and extended fits.

Table 1 Station-specific coefficients K_station obtained from Eq. (1) and effective coefficients K_eff(h)=Kn(h)^a obtained from the joint density-dependent fit.

For the combined dataset, the coefficient of determination increases from R²≈0.90 for the fit of Eq. (1) to R²≈0.99 for the density-dependent formulation. Figure 2 compares the effective avalanche lengths obtained from the CORSIKA simulations with the corresponding model predictions. The density-dependent formulation shows better agreement with the simulations, particularly at larger avalanche lengths where the fit of Eq. (1) shows larger deviations from the 1:1 relation.

Figure 3 Electron (a) and gamma-ray (b) energy spectra for detector levels located 25, 50, 100, and 200 m below the lower boundary of the electric field for the Aragats configuration. The spectra are normalized per injected seed electron.

Figure 4 Electron (a) and gamma-ray (b) counts at detector levels located 25, 50, 100, and 200 m below the lower boundary of the electric field for the Aragats configuration. Particle counts are normalized per injected seed electron.

4 Particle Propagation beyond the AEF

To evaluate the dependence of particle energy spectra and total particle counts on the detector distance, additional simulations were performed for the Aragats configuration. All simulation parameters were kept unchanged, and only the detector position was varied. The detector was placed 25, 50, 100, and 200 m below the lower boundary of the electric field.

The resulting electron and gamma-ray energy spectra are presented in Fig. 3.

As the field-free propagation distance increases, the electron spectrum is progressively attenuated, with a pronounced suppression of the high-energy component. The gamma-ray component persists over larger distances, although its intensity also decreases with increasing propagation distance. This different transport behavior becomes particularly evident when the corresponding particle counts are considered (Fig. 4).

The electron population decreases rapidly and monotonically with increasing field-free propagation distance, whereas the gamma-ray population is attenuated more gradually. The extension of the simulations to 200 m confirms that the gamma-ray component also undergoes substantial attenuation at larger propagation distances. A separation of 100 m represents a configuration in which the electron component is already strongly attenuated while a substantial gamma-ray component is still preserved.

For the following analysis, E^max is defined as the highest energy for which at least three consecutive energy bins contain more than five counts. This criterion reduces sensitivity to isolated high-energy outliers. The electron and gamma-ray energies are denoted by $E_{\mathrm{e}}^{\max }$ and $E_{\mathit{\gamma }}^{\max }$, respectively. The approximate gap length (free path distance) can then be estimated as:

$$\begin{array}{}\text{(11)}& \mathrm{FPD}={\displaystyle \frac{E_{\mathrm{e}}^{\max }\left(\mathrm{exit}\right)-E_{\mathrm{e}}^{\max }\left(\det \right)}{C_{\mathrm{2}}\left(h\right)}},\end{array}$$

where C₂(h) is the mean energy loss of the electron at altitude h. However, $E_{\mathrm{e}}^{\max }$(exit) is not directly available from detector measurements. Instead, the maximum gamma-ray energy $E_{\mathit{\gamma }}^{\max }$(det) measured at the detector provides an observable quantity that can be related to the particle energies at the field boundary. Unlike electrons, which undergo continuous energy losses during propagation through air, gamma rays propagate through stochastic interactions whose probabilities depend on photon energy, including Compton scattering, photoelectric absorption, and electron–positron pair production. To quantitatively describe this propagation beyond the accelerating electric-field region, Chilingarian et al. (2021) introduced an empirical relation for estimating the free path distance (FPD). Since $E_{\mathrm{e}}^{\max }$(exit) is not directly available from detector measurements, it is estimated using the empirical relation $E_{\mathrm{e}}^{\max }\left(\mathrm{exit}\right)\approx C_{\mathrm{1}}{E}_{\mathit{\gamma }}^{\max }\left(\mathrm{exit}\right)$, where C₁ relates the characteristic upper energies of the electron and gamma-ray components at the field boundary. Substitution of this relation into Eq. (11) gives

$$\begin{array}{}\text{(12)}& \mathrm{FPD}={\displaystyle \frac{C_{\mathrm{1}}{E}_{\mathit{\gamma }}^{\max }\left(\mathrm{exit}\right)-E_{\mathrm{e}}^{\max }\left(\det \right)}{C_{\mathrm{2}}\left(h\right)}}.\end{array}$$

Here $E_{\mathit{\gamma }}^{\max }$(exit) represents the characteristic gamma-ray energy at the exit of the electric-field region. For the 100 m separation used to calibrate the FPD relation, $E_{\mathit{\gamma }}^{\max }\left(\mathrm{exit}\right)\approx E_{\mathit{\gamma }}^{\max }\left(\det \right)$, as demonstrated by the simulated spectra below. $E_{\mathrm{e}}^{\max }$(det) is the characteristic electron energy at the detector level, C₁ is an empirical electron–gamma energy coefficient, and C₂(h) is the altitude-dependent electron energy loss rate in air. The coefficient C₁ is determined jointly for the four stations. For this analysis, the field-to-detector separation was fixed at 100 m for all four stations. The electric field strengths were selected for each station as follows: Aragats (210 kV m⁻¹), LHAASO (190 kV m⁻¹), Lomnický Štít (220 kV m⁻¹), and Nor Amberd (240 kV m⁻¹). Electron and gamma-ray spectra were obtained at both the exit and detector levels. A characteristic maximum energy E^max was defined as the highest energy for which at least three consecutive energy bins contain more than five counts. This criterion reduces sensitivity to high-energy outliers.

Gamma-ray spectra change little over the 100 m propagation distance, with characteristic energies of $E_{\mathit{\gamma }}^{\max }\left(\mathrm{exit}\right)\approx \mathrm{29}$–30 MeV and $E_{\mathit{\gamma }}^{\max }\left(\det \right)\approx \mathrm{30}$–31 MeV, supporting the approximation $E_{\mathit{\gamma }}^{\max }\left(\mathrm{exit}\right)\approx E_{\mathit{\gamma }}^{\max }\left(\det \right)$. Accordingly, the detector value is used as a proxy for the corresponding field-exit quantity in Eq. (12). Conversely, electron spectra show a clear shift toward lower energies. At Aragats, the electron maximum energy dropped from approximately 38 MeV at exit to around 16 MeV at the detector after 100 m. This behavior is shown in Figs. 5 and 6.

Figure 5 The energy spectra of electrons at the exit and detector levels for all four stations: (a) Aragats (210 kV m⁻¹), (b) LHAASO (190 kV m⁻¹), (c) Lomnický Štít (220 kV m⁻¹), and (d) Nor Amberd (240 kV m⁻¹). The vertical separation between exit and detector is fixed at 100 m in all cases.

Figure 6 Gamma-ray energy spectra at the exit and detector levels for all four stations: (a) Aragats (210 kV m⁻¹), (b) LHAASO (190 kV m⁻¹), (c) Lomnický Štít (220 kV m⁻¹), and (d) Nor Amberd (240 kV m⁻¹). The vertical separation between the exit and detector levels is fixed at 100 m in all cases. The characteristic high-energy endpoints show little change over the propagation distance, and the effective maximum energies are indicated by vertical dashed lines.

The relation between the characteristic electron and gamma-ray energies is parameterized by the coefficient C₁ in Eq. (12). The proportionality coefficient C₁ is derived directly from the exit and detector spectra using a bootstrap-based statistical procedure. For each station, the exit energy histograms of gamma rays and electrons were taken as the starting point. A Poisson bootstrap procedure was then applied, in which each energy bin count N_i was resampled according to

$$\begin{array}{}\text{(13)}& N_i^{\left(\mathrm{b}\right)}\left(E\right)\sim \mathrm{Poisson}\left(N_i\right),\end{array}$$

where b denotes the bootstrap realization. A total of B=2000 bootstrap realizations were used. For each bootstrap realization, the characteristic energies E_e(exit)^(b), E_e(det)^(b) and E_γ(exit)^(b) were recalculated using the same three-consecutive-bin criterion defined above. Then, the corresponding coefficient was calculated as

$$\begin{array}{}\text{(14)}& C_{\mathrm{1}}^{\left(\mathrm{b}\right)}=\left(\mathrm{100}\cdot C_{\mathrm{2}}\left(h\right)+E_{\mathrm{e}}^{\left(\mathrm{b}\right)}\left(\det \right)\right)/E_{\mathit{\gamma }}^{\left(\mathrm{b}\right)}\left(\mathrm{exit}\right).\end{array}$$

The bootstrap realizations yield a C₁ distribution for each station. The median values obtained for individual stations are:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{Aragats}:C_{\mathrm{1}}\approx \mathrm{1.31},\\ {\displaystyle }& {\displaystyle }\mathrm{LHAASO}:C_{\mathrm{1}}\approx \mathrm{1.39},\\ {\displaystyle }& {\displaystyle }\mathrm{Lomnick}\mathrm{ý}\mathrm{Š}\mathrm{t}\mathrm{í}\mathrm{t}:C_{\mathrm{1}}\approx \mathrm{1.46},\\ {\displaystyle }& {\displaystyle }\mathrm{Nor}\mathrm{Amberd}:C_{\mathrm{1}}\approx \mathrm{1.44}.\end{array}$$

In the present formulation, C₁ is calibrated for a field-to-detector separation of 100 m. The station-specific C₁ distributions show typical spreads of approximately ±0.03–0.05. Combining the bootstrap realizations from all four stations yields a pooled median of C₁≈1.38, while a global least-squares optimization, minimizing the deviations of the predicted FPD from the 100 m propagation scale, gives C₁≈1.375. The close agreement between these independent estimates supports the stability of the inferred coefficient across the atmospheric conditions considered here. Accordingly, C₁≈1.37 is adopted for the subsequent FPD parameterization.

The second parameter of the model, C₂(h), characterizes the average energy loss of electrons per unit path length, which depends explicitly on atmospheric density. For each station, a representative energy for calculation losses was defined as

$$\begin{array}{}\text{(15)}& E^{*}=\left(E_{\mathrm{e}}\left(\mathrm{exit}\right)+E_{\mathrm{e}}\left(\det \right)\right)/\mathrm{2},\end{array}$$

which approximates the typical electron energy during propagation. The corresponding mass stopping power was obtained from the NIST ESTAR database (Berger et al., 2005) using linear interpolation and multiplied by the local air density to obtain the linear stopping power.

The resulting values are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{Aragats}:C_{\mathrm{2}}\approx \mathrm{0.20}\mathrm{MeV}{\mathrm{m}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }\mathrm{LHAASO}:C_{\mathrm{2}}\approx \mathrm{0.17}\mathrm{MeV}{\mathrm{m}}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }\mathrm{Lomnick}\mathrm{ý}\mathrm{Š}\mathrm{t}\mathrm{í}\mathrm{t}:C_{\mathrm{2}}\approx \mathrm{0.21}\mathrm{MeV}{\mathrm{m}}^{-\mathrm{1}},\end{array}$$

$$\mathrm{Nor}\mathrm{Amberd}:C_{\mathrm{2}}\approx \mathrm{0.23}\mathrm{MeV}{\mathrm{m}}^{-\mathrm{1}}.$$

The FPD was then calculated according to Eq. (12), where C₁ was derived from spectral data, C₂(h) from stopping-power tables, and the energies from the simulated spectra.

Applying this expression to the simulation data for four stations yields FPD values ranging from 85 to 105 m, in agreement with the expected ∼100 m scale.

To compare the empirical coefficient C₁=1.2 with the optimized coefficient C₁≈1.37, FPD distributions were calculated from the bootstrap realizations for all four stations using Eq. (12). The resulting distributions are shown in Fig. 7, with the mean values and corresponding 68 % intervals indicated.

Figure 7 Comparison of free path distance (FPD) distributions obtained using the empirical coefficient C₁=1.2 and the optimized value C₁=1.37 for four high-altitude stations (Aragats, LHAASO, Lomnický Štít, Nor Amberd). The distributions are derived from bootstrap resampling of the simulated spectra. Colored violins show the density of FPD values, while black markers denote mean values with 16th–84th percentile intervals. The dashed line indicates the reference field-to-detector separation of 100 m.

As shown in Fig. 7, the empirical coefficient C₁=1.2 systematically underestimates the FPD at all four stations, whereas the optimized value C₁≈1.37 gives values closer to the 100 m propagation distance. The remaining spread reflects the statistical variability of the high-energy spectral tails.

To test the robustness of the FPD estimate, 50 additional independent CORSIKA simulations were performed for the Aragats station at an electric field of 210 kV m⁻¹, with a 50 m distance between the lower boundary of the electric-field region and the detector level. The resulting data were used to recalculate the FPD using the same coefficients obtained for the 100 m separation adopted for calibration of the FPD relation. For the 50 m distance, the median of the reconstructed FPD distribution was 69.9 m, with a 16th–84th percentile interval of 58.65–77.35 m:

$$\begin{array}{}\text{(16)}& \mathrm{FPD}={\mathrm{69.9}}_{-\mathrm{11.2}}^{+\mathrm{7.5}}\mathrm{m}.\end{array}$$

The reconstructed FPD is shifted toward larger values by approximately 20 m. Nevertheless, the reconstructed distribution correctly identifies the short-distance range over which electrons emerging from the electric-field region can reach the detector level.

The obtained 16th–84th percentile interval should be regarded as a lower bound on the expected uncertainty, since the additional uncertainty associated with the experimental procedure of recovering energy spectra from energy-release histograms is not included in the simulations.

5 Conclusions

In this work, CORSIKA simulations were used to investigate two complementary stages of relativistic runaway electron avalanche development: avalanche growth within the atmospheric electric field and subsequent particle propagation beyond the accelerating region. The simulations were performed for four high-altitude observational sites spanning different atmospheric densities and electric-field conditions.

The effective avalanche lengths derived from the vertical electron profiles show station-to-station differences when a single coefficient is fitted to the combined dataset. The station-specific coefficients range from $K=(\mathrm{8.6}\pm \mathrm{0.2})\times {\mathrm{10}}^{\mathrm{3}}$ to $(\mathrm{11.2}\pm \mathrm{0.3})\times {\mathrm{10}}^{\mathrm{3}}$ kV. Introducing an additional density-dependent term reduces the inter-station scatter, increasing the coefficient of determination from R²≈0.90 to R²≈0.99. The joint fit yields $K\approx (\mathrm{6.39}\pm \mathrm{0.26})\times {\mathrm{10}}^{\mathrm{3}}$ kV and $a=-\mathrm{0.923}\pm \mathrm{0.083}$, indicating an additional approximately inverse dependence on atmospheric density beyond that already included in the runaway-threshold scaling.

The simulations also characterize the different evolution of the electron and gamma-ray components after leaving the accelerating electric-field region. Electron spectra and particle counts decrease rapidly with field-free propagation distance, whereas the gamma-ray component is attenuated more gradually. Within the FPD formulation, the characteristic electron and gamma-ray energies at the field boundary are related through the coefficient C₁. Calibration for a 100 m field-to-detector separation gives C₁≈1.37. The density-dependent electron energy-loss coefficient C₂(h), obtained from ESTAR stopping-power data, ranges from approximately 0.17 to 0.23 MeV m⁻¹ across the four stations. Using these parameters, the reconstructed FPD values range from approximately 85 to 105 m, consistent with the 100 m propagation scale used for the calibration.

These results are derived for the idealized uniform electric-field configurations considered here and should not be generalized to arbitrary thunderstorm electric-field structures without further validation.