Threshold atmospheric electric fields for initiating relativistic runaway electron avalanches: theoretical estimates and CORSIKA simulations

We examine the threshold atmospheric electric field (E_th) required to initiate a runawayavalanche in Earth's atmosphere. We compare the traditional, thirty-year-old theoretical threshold with its recently updated value and with the threshold derived from CORSIKA-simulated avalanches (E_z). The altitude dependence of these thresholds is analyzed, considering how changes in air density affect avalanche development. This study is vital for understanding high-energy atmospheric phenomena in both the lower and upper atmosphere, including thunderstorm ground enhancements (TGEs) and gamma glows, and for refining atmospheric electric field (AEF) models based on particle flux measurements.

1 Introduction

Free electrons are abundant in the troposphere. The altitude at which their flux reaches its highest point, called the Regener–Pfotzer maximum (Regener, 1933). It depends on the geomagnetic cutoff rigidity (Rc), the type of particles being measured, and the phase and strength of the solar cycle. Recent observations, supported by PARMA4.0 calculations (Sato, 2016), show that at middle to low latitudes (Rc = 3–8 GV), the highest flux of charged particles occurs at altitudes around 12–14 km (see Fig. 3 in Ambrozová et al., 2023).

Atmospheric electric fields (AEFs) generated by thunderstorms transfer energy to free electrons, accelerate them, and, under certain conditions, induce electron-photon avalanches. Gurevich et al. (1992) identified the conditions necessary for extensive multiplication of electrons from an energetic seed electron injected into a strong AEF region. This process is known as the Relativistic Runaway Electron Avalanche (RREA; Babich et al., 2001; Alexeenko et al., 2002). A numerical approach for solving the relativistic Boltzmann equation for runaway electron beams (Symbalisty et al., 1998) aids in estimating the threshold AEF (Babich et al., 2001; Dwyer, 2003) required to trigger RREA. As demonstrated by GEANT4 and CORSIKA simulations (Chilingarian et al., 2012, 2022, 2024), the RREA process is a threshold phenomenon, with avalanches initiating when the atmospheric AEF exceeds a certain threshold that depends on air density. The AEF must also be sufficiently extended to support the growth of avalanches. At standard temperature and pressure in dry air at sea level, E_th≈ 2.80 kV cm${}^{-\mathrm{1}}\times n$, where air density n is relative to the International Standard Atmosphere (ISA) sea-level value (see the recent update of the threshold energy E_th≈ 2.67 kV cm${}^{-\mathrm{1}}\times n$ in Dwyer and Rassoul, 2024). This threshold field is slightly higher than the breakeven field, which corresponds to theelectron energy at which minimum ionization occurs. If electrons traveled exactly along AEF lines, it would define the threshold for runaway electron propagation and the start of avalanche formation. However, the paths of electrons deviate due to Coulomb scattering with atomic nuclei and Møller scattering with atomic electrons, causing deviations from the near-vertical AEF. Additionally, secondary electrons produced by Møller scattering are not generated along the field line; therefore, AEFs must exceed the theoretical RREA threshold E_th by approximately 10 %–20 % for electrons to run away and trigger an avalanche.

2 CORSIKA simulations of RREAs reaching the high-altitude stations

To understand how avalanches develop in an electrified atmosphere and to compare the new and updated E_th with the particle-intensity abrupt growth, we used the CORSIKA code (Heck et al., 1998), version 7.7500, which accounts for the effect of AEFs on particle transport (Buitink et al., 2009). The growth of RREA increases the cloud's electrical conductivity. Numerous studies (Marshall et al., 1995; Stolzenburg et al., 2007) have indicated that lightning flashes tend to occur when the applied electric field exceeds the RREA threshold by roughly 20 %–30 %. RREA simulation codes do not include a lightning initiation mechanism. Therefore, one can artificially raise the AEF strength beyond a realistic value to produce billions of avalanche particles; however, this approach lacks physical justification. As a result, we do not test AEFs stronger than 2.5 kV cm⁻¹ at altitudes of 3–6 km. The RREA simulation was performed for vertical seed electrons with a uniform AEF that exceeded the E_th by a few tens of percent. An introduced fixed uniform AEF shifts the surplus to E_th at different heights by different percentages, corresponding to air density. The chosen seed electron energy spectrum was based on the EXPACS Excel-based program (Sato, 2016), following a power-law with an index of 1.173 for energies from 1 to 300 MeV. During TGE events on Aragats, the typical distance to the cloud base is estimated to be 25–200 m (see Fig. 17 in Chilingarian et al., 2020); therefore, in our simulations, particle propagation continued in dense air for an additional 25, 50, 100, and 200 m before detection.

The simulations included 1000 to 10 000 events for AEF strengths from 1.55 to 2.5 kV cm⁻¹. Electron and gamma-ray propagation were tracked until their energies dropped to 0.05 MeV. Each simulation event corresponds to the propagation of a single seed electron; multiple events were used to obtain statistically stable averages to reliably estimate the resulting threshold electric fields. The CORSIKA code models the development of RREA by calculating the number of electrons and gamma rays at different stages of the cascade development at 200-meter intervals. At all stations, the atmospheric electric field was implemented as a vertically uniform layer with a thickness of 2000 m above the observation levels. Besides the Aragats and Nor Amberd research stations on the slopes of Mt. Aragats in Armenia, we also conducted simulations for the Slovakian and Chinese research stations at Lomnický Štít (Chum et al., 2020) and the Tibetan plateau. LHAASO (Large High Altitude Air Shower Observatory, Aharonian et al., 2023) is situated at 4410 m a.s.l. It provides an ideal platform for studying atmospheric particle acceleration, owing to its thin atmosphere and the high likelihood of runaway electron avalanche formation. For LHAASO, we present CORSIKA simulation results showing increases in electron and photon fluxes under AEF strengths ranging from 1.55 to 1.9 kV cm⁻¹. The number of electrons and photons was recorded at altitudes ranging from 6510 to 4510 m. Lomnický Štít is located at an altitude of 2630 m in Slovakia. CORSIKA simulations were performed for various vertical AEFs ranging from 1.9 to 2.3 kV cm⁻¹. The number of electrons and photons was recorded at altitudes ranging from 4734 to 2734 m. Significant increases in flux were observed with stronger-than-threshold fields, confirming the development of robust RREA. Saturation trends in the growth of electrons and photons suggest that the threshold field, E_th, at Lomnický Štít is approximately 2.3 kV cm⁻¹. These results support earlier findings from Aragats and Nor Amberd and emphasize the altitude dependence of E_th. Due to the thinner air at LHAASO, the TGEs occurred at a much lower value of 1.7 kV cm⁻¹. In Fig. 1a–d, we display the development of RREA at different atmospheric depths and for various physically justified strengths of the AEF. The curves are scaled for a single seed electron for easier comparison with experimentally measured intensities.

For large values of AEFs, the number of avalanche particles rose exponentially. For lower values of AEF, we observe saturation of the particle flux when AEF becomes lower than the threshold electric field (dependent on air density); the RREA process attenuates before reaching the observation level (see the yellow and blue curves in Fig. 1a–d).

Figure 1 Longitudinal development of relativistic runaway electron avalanches (RREA) at four high-altitude observation sites: (a) Aragats, (b) Lomnický Štít, (c) Nor Amberd, and (d) LHAASO. The number of electrons is normalized to a single seed electron and shown as a function of depth within the electric field region. For each site, simulations were performed for several electric field strengths, as indicated in the legends. Avalanche development is sampled every 200 m across the field region. After exiting the electric field, electron propagation is followed for an additional 100 m in free space before reaching the detector.

Figure 1 illustrates the dependence of electron multiplication on electric field strength and highlights the altitude-dependent conditions required for sustained avalanche development.

We estimate the “simulated” thresholds, E_z values, at the heights where the number of avalanche particles stops rising, as shown in Fig. 2.

Figure 2 Simulated threshold electric field strength, E_z, versus altitude for several high-altitude stations. The threshold is defined as the electric field strength at which the growth of avalanche electrons saturates.

Figure 3 Simulated threshold electric field strength, E_z, and theoretical threshold electric fields, E_th, as a function of saturation (rise stopping) altitude.

Table 1 Excess of E_z over E_th. Stopping altitudes and theoretical threshold field comparisons for heights 2500–5550 m.

In Fig. 3 and Table 1, we compare the “simulated” threshold E_z with the theoretical ones. Simulations derive higher values than theoretical estimates, especially for high E_th values (low altitudes) at all four research stations. Theoretical threshold fields are computed as 2.67 and 2.80 kV cm${}^{-\mathrm{1}}\times n$. The percentage of enhancement indicates how much the applied field exceeds the theoretical thresholds. Strong AEFs, where the cascade did not attenuate, were not included in the table.

3 Discussion and conclusion

Both the classical threshold field (E_th≈ 2.80 kV cm${}^{-\mathrm{1}}\times n$) and its updated version (E_th≈ 2.67 kV cm${}^{-\mathrm{1}}\times n$) are derived under idealized assumptions; the difference between them results from refinements in modeling particle energy los sprocesses. The earlier estimate of 2.80 kV cm${}^{-\mathrm{1}}\times n$ was based on basic energy balance considerations using older ionization loss models and assumed monoenergetic electrons. This threshold is slightly above the breakeven field, where energy gain equals average energy loss. The updated 2.67 kV cm${}^{-\mathrm{1}}\times n$ value, introduced by Dwyer and Rassoul (2024), incorporates more accurate relativistic Boltzmann solutions, improved ionization and bremsstrahlung cross-sections, and a probabilistic treatment of runaway thresholds across realistic energy spectra. While both thresholds assume idealized, field-aligned electron motion in a uniform medium, the updated value is physically more consistent. It predicts a slightly lower field strength needed for initial runaway. However, CORSIKA simulations show that this refined threshold is insufficient for sustained avalanche growth under real atmospheric conditions due to scattering and finite path effects. Moreover, it deviates more from the simulated value than the “classical”, 30-year-old estimate. Multiple physical processes inhibit ideal runaway propagation. Coulomb scattering with atmospheric nuclei and Møller scattering with electrons cause substantial angular deflection and energy redistribution. Secondary electrons are not generated strictly along the field direction, and many lose energy before gaining sufficient momentum to continue avalanche growth. As a result, electrons must be accelerated in fields stronger than the threshold to overcome these losses and maintain avalanche conditions. CORSIKA simulations, which incorporate all major interaction mechanisms – includingCoulomb and Møller scattering, bremsstrahlung losses, finite propagation distances, andrealistic secondary cosmic ray spectra – show that avalanches fully developonly when the applied field exceeds the theoretical threshold by a measurable margin. For the updated 2.67 kV cm${}^{-\mathrm{1}}\times n$ value, we observe a required excess of approximately 20 %–22 % at the Aragats station (∼ 3200–4200 m a.s.l.), whereas for the classical 2.80 kV cm${}^{-\mathrm{1}}\times n$ threshold, the excess is typically 15 %–17 %. Interestingly, this required excess decreases with increasing air density, as observed inthe Nor Amberd simulations. At lower altitudes (∼ 2500–2700 m a.s.l.), the differencebetween the applied and threshold fields is reduced: only 14 %–16 % above 2.67 kV cm${}^{-\mathrm{1}}\times n$, and about 9 %–11 % above 2.80 kV cm${}^{-\mathrm{1}}\times n$.

This trend can be explained as follows: In denser air, the chances of energy-loss interactions increase, but so does the likelihood of electron multiplication through ionization and bremsstrahlung over shorter distances. The avalanche can develop more quickly because seed electrons encounter more target atoms in a given path length. As a result, the necessary “headroom” above the threshold field for sustained multiplication is smaller. Simply put, the efficiency of avalanche formation improves in denser air, even though the absolute threshold field is higher. This results in a smaller relative excess being required above the theoretical threshold. Therefore, although the threshold field scales linearly with air density, the required enhancement factor does not. It decreases with increasing density due to a balance between energy loss and multiplication processes, all of which are faithfully captured in the CORSIKA simulation framework. This emphasizes the importance of altitude-dependent analysis in interpreting Thunderstorm Ground Enhancements (TGEs) and suggests that scaling laws based solely on density may overlook subtler effects arising from atmospheric structure and shower development dynamics.Among these effects is the local temperature profile, which can modify air density and, consequently, slightly affect the effective threshold field. A more detailed treatment incorporating measured or modeled temperature profiles could further refine threshold estimates for individual events; however, such event-specific modeling is beyond the scope of the present work.