Calculation of beams of positrons, neutrons, and protons associated with terrestrial gamma ray flashes

2.1.1 Stepped Lightning Leader

For the stepped lightning leader, we use the model of Xu et al. [2012a] with the same parameters. The leader is vertical, 4 km long, and has a tip radius of 1 cm. It is assumed to be equipotential and embedded in an external thundercloud electric field of 0.5 kV/cm, far from the leader. This model approximates the leader at the moment of stepping: the space stem has connected to the leader, and the full electric potential is now on the new leader tip. At this moment the leader tip “explodes” with ionization, similarly as described by Sun et al. [2013], creating an inception cloud [Briels et al., 2008] and later a streamer corona. The field of the leader is tested by inserting 50 electrons with 0.1 eV energy on the symmetry axis 30 cm ahead of the leader tip, as in the work of Xu et al. We also follow the approach of Xu et al. [2012a] by not taking the space charge effects of the developing corona discharge into account, but only that of the stationary leader; this approximation will be justified in Figure 3 for the electrons with the highest energy. Therefore, the approximation of taking only the electrostatic leader field into account is reasonable, even if the inception cloud develops into a relativistic impact ionization front [Luque, 2014].

Rather than approximating the leader as a cylinder with semispherical caps as Xu et al. [2012a], we approximate it as an ellipsoid with a length of 4 km and a curvature radius of 1 cm at the tip. This has the advantage that we can calculate the electric field E(r) analytically, as summarized in Appendix A. Figure 1 shows the electric field strength in the vicinity of the leader tip. It shows that the ellipsoid is a reasonable approximation when comparing with Figure 1a of Xu et al. [2012a], and it illustrates the strong field enhancement close to the leader and its tip. The field is approximately 500 kV/cm at 30 cm ahead of the leader tip, thus certainly large enough to accelerate the electrons into the run-away regime.

2.1.2 Air Composition and Monte Carlo Approach

We model the air as consisting of 78.12% N₂, 20.95% O₂, and 0.93% Ar. To control the air density as a function of altitude, we use the barometric formula with a scale height of 8.33 km. We assume the upper leader tip where the electrons are accelerated to lie at 16 km altitude which corresponds to an air density of 1/10 of the density at sea level if the temperature change with altitude is taken into account.

We trace the positions of electrons and photons in three dimensions with a Monte Carlo code where the neutral air molecules are treated as a random background with appropriate statistical weight. Between collisions, the electrons follow classical or relativistic trajectories within the given electric field, depending on their energy, while the photons move with the speed of light into the direction of emission. The collisions are treated with a Monte Carlo scheme.

We now list the collision types included.

2.1.3 Cross Sections for Electrons

Details on our Monte Carlo code and its validation can be found in Köhn et al. [2014] and Köhn and Ebert [2014b]. There we also describe which collisions we take into account and how we have implemented them into our code. We remark here that the choice of cross sections is essential for correct results as shown in Appendix C of Köhn and Ebert [2014b] and that there is continued research on the cross sections for the propagation and production of high-energy particles in the atmosphere. Of particular importance is the correct modeling of Bremsstrahlung that creates the X and gamma rays from the high-energy electrons.

2.1.4 Electron-Nucleus Bremsstrahlung

Electron-nucleus Bremsstrahlung is the process when a free electron scatters at a nucleus and emits a photon, electron-electron Bremsstrahlung is the same process where the electron scatters at another electron and emits a photon. As electron-electron Bremsstrahlung is frequently considered as negligible compared to electron-nucleus Bremsstrahlung, the general term Bremsstrahlung refers typically to electron-nucleus Bremsstrahlung.

Different cross-sections for electron-nucleus Bremsstrahlung are used in different databases and by different researchers. In the field of terrestrial gamma ray flashes, Carlson et al. [2010] use the Geant4 simulation tool kit with its intrinsic cross sections [Agostinelli et al., 2003]. Dwyer [2007] uses the Bethe-Heitler cross section [Bethe and Heitler, 1934; Heitler, 1944] resolving the full geometry of the process; he includes an atomic form factor to take the structure of the atomic shell into account. Xu et al. [2012a] use the simple product ansatz of Lehtinen [2000] to relate the energy and the direction of the emitted photons.

Now Geant4 [Agostinelli et al., 2003] uses cross sections appropriate for the large atomic numbers Z of heavy nuclei while nitrogen and oxygen have Z = 7 and 8; thus, the cross sections implemented in Geant4 are not appropriate to simulate the production of Bremsstrahlung photons in air. According to Shaffer et al. [1996], the old Bethe-Heitler theory keeps being the appropriate theory for Z <29 for electron energies between 1 keV and 1 GeV, hence for the production of TGFs, as we have already discussed in Köhn and Ebert [2014a]. In Appendix F of Köhn and Ebert [2014a] we have also shown that the atomic form factor used by Dwyer [2007] is close to unity in the relevant cases and thus negligible. The product ansatz of Lehtinen [2000] that is used by Xu et al. [2012a] is not compatible with the full quantum field theoretical treatment of collisions where the photon obtains almost all the electron energy, as we have discussed in Appendix E of Köhn and Ebert [2014a]; we will compare results of TGF calculations under the same conditions using either the cross sections of Lehtinen [2000] or of Köhn and Ebert [2014a] in 13.

We use the doubly differential cross section derived by Köhn and Ebert [2014a] for the relation between the photon energy and the angle Θ_i between the incident electron and the emitted photon. This cross section has been obtained by integrating the direction of the emitted electron out in the Bethe-Heitler cross sections, and it is implemented using rejection sampling as described in Knuth [1997]. The scattered electron keeps its initial direction.

2.1.5 Electron-Electron Bremsstrahlung

As databases like Geant 4 concentrate on the Bremsstrahlung for metals like iron (Z = 26) or lead (Z = 82), electron-electron Bremsstrahlung is considered as irrelevant, because it scales with Z rather than with Z². Furthermore, the photons emitted in electron-electron Bremsstrahlung by nitrogen or oxygen (Z = 7,8) are negligible as well compared to electron-nucleus Bremsstrahlung. However, we have shown recently [Köhn et al., 2014] that electron-electron Bremsstrahlung also ejects the shell electron on which the free electron is scattered and that these electrons have higher energies than those created by normal impact ionization [Yong-Ki and Santos, 2000]. Using the cross sections of Tessier and Kawrakow [2007], the effect is important in air for electron energies of several MeV. Electron-electron Bremsstrahlung hence largely increases the number of electrons with energies above 1 MeV, and it subsequently contributes substantially to the number of high-energy photons from a negative stepped lightning leader. We remark here that electron-electron Bremsstrahlung also has been included recently in simulation packets like CORSIKA [2013] and EGS5 [2005] simulating extensive air showers, as well as in EGSnrc [2014] for medical applications, but not yet in Geant4 [Agostinelli et al., 2003] to the best of our knowledge.

2.1.6 Cross Sections for Photons

Figure 2 shows cross sections of photon processes as a function of the photon energy. For the photons, we use the cross sections for photoionization (Evaluated Photon Interaction Data (EPDL) Database, Photon and Electron Interaction Data, 1997, http://www-nds.iaea.org/epdl97/), Compton scattering [Greiner and Reinhardt, 1995; Peskin and Schroeder, 1995], hadron production [Fuller, 1985], and pair production (Bethe-Heitler in the integrated form of Köhn and Ebert [2014a]). Figure 2 shows that photoionization is dominant for photon energies below 1 keV. Thus, new electrons are created not only by electron impact ionization and the electron-electron Bremsstrahlung process, but also by photoionization.

2.2.1 Distribution of Electrons in Energy and Space

In Köhn et al. [2014], we have already presented the electron energy distribution ahead of the stepped lightning leader after 24 ns of evolution when either including or neglecting the electron-electron Bremsstrahlung. There, we have presented the electron energy distribution for different runs with different realizations of random numbers, and we have seen that the energy distribution is stable against different sets of random numbers already for 50 initial particles; hence, 50 initial electrons is already enough for good statistics. Due to the limited space of a Fast Track Communication, we could not present the buildup of the spatial distribution. Therefore, we present here in Figure 3 this evolution including electron-electron Bremsstrahlung when 50 test electrons are inserted 30 cm ahead of the leader tip. The leader is indicated in black. The spatial distributions of the electrons after 5 ns, 10 ns, 15 ns, and 20 ns are plotted; the continuous colors indicate the electron densities with energies below 1 MeV, the color lines the electron densities with energies above 1 MeV in the xz plane where a slice was evaluated in the y direction from −3 cm to 3 cm.

The figure shows that at all instances, the high-energy electrons are ahead of the lower energy electrons and more on axis. The reason is obviously that the high-energy electrons are accelerated continuously while the lower energy electrons have lost energy in collisions and these collisions also lead to a widening of the electron beam. The figure shows as well that new ionization patches are created at the sides of the leader, probably due to photoionization created by Bremsstrahlung photons. It should be noted though that the motion of the low-energy electrons is not quite physical as we neglect the space charge effects of the newly created ionization, just like Xu et al. [2012a, 2012b]. However, as there is a clear spatial separation between the electron populations at different energies, we argue that the calculation approximates the high-energy spectrum of the electrons well. We remark that electrons with energies of 10 keV, 100 keV, 1 MeV, and 10 MeV move with 19.5, 54.8, 94.1, and 99.9% of the speed of light and that light travels 6 m within 20 ns, while the fastest electrons are ≈ 4 m from the point of injection after 20 ns; see Figure 3d. The electrons with energies above 1 MeV are concentrated in one region on axis after 5 and 10 ns; after 15 ns, new patches with high-energy electrons have formed in new beam directions slightly off axis. We emphasize that the channel-like structures forming from 10 ns on are not streamers, as space charge effects are not included, but they are probably rather ionization traces of the created high-energy electrons, similarly as in cosmic particle showers [Köhn and Ebert, 2014b], but enhanced by the electric field.

We remark that the geomagnetic field has no influence on the electrons at these altitudes, and we will discuss the role of the geomagnetic field in more detail in 16.

2.2.2 Distribution of Photons in Energy and Space

Figure 4 shows the photon energy distribution after 24 ns. As we have found in Köhn et al. [2014], the spectrum of Bremsstrahlung photons, and in particular the high energy tail, does not change considerably after about 15 ns. For reasons of computer memory, we anyhow have inserted only 50 test electrons into the leader field and we follow the motion of the electrons until 24 ns. The total number of photons with energies between 0.01 eV and 10 MeV in our simulation is approximately 5000; thus, there are 100 photons per initial electron. As we have neglected the space charge effects of the developing corona discharge, the low-energy part of the distribution is not quite physical and we concentrate here on photon energies above 10 keV. The maximal photon energy after 24 ns is approximately 10 MeV. Figure 4 also shows that for energies above 1 MeV, the distribution can be fitted well by the exponential .

From our analysis in Köhn and Ebert [2014a], we know that photons with energies above 1 MeV are emitted predominantly in forward direction relative to the direction of the incident electron. Since the electrons, which can create such photons, mainly move upward (see Figure 3), this suggests to assume that the beam of initial photons is monodirectional. We will use this assumption in 22 as a test case to simulate the production of positrons and hadrons.

2.2.3 Results for Different Bremsstrahlung Cross Sections

We tested the dependence of the simulation results on different Bremsstrahlung cross sections. Figure 5 shows the photon energy distribution after 24 ns for energies above 1 keV. The two curves are both without electron-electron Bremsstrahlung, and either the product ansatz of Lehtinen [2000] or the integrated Bethe-Heitler cross section of Köhn and Ebert [2014a] for the electron-nucleus Bremsstrahlung is used. The plot shows that the product ansatz of Lehtinen [2000] that is used by Xu et al. [2012a, 2012b] substantially overestimates the number of photons with energies above 20 keV, hence also the photons in the MeV range. We note here that there are some fluctuations due to the small number of produced photons; still the tendency of overestimating the number of high-energy photons is clearly observed.

3.1.1 Compton Scattering

The cross section dσ/dΩ for Compton scattering is [Greiner and Reinhardt, 1995; Peskin and Schroeder, 1995]

(1)

where α ≈ 1/137 is the fine structure constant, h ≈ 6.6 · 10⁻³⁴ Planck's constant, m ≈ 9.1 · 10⁻³¹ kg the electron mass, c ≈ 3 · 10⁸ m/s the speed of light, and

(2)

is the energy of the scattered photon as a function of the incident photon energy E_γ and of the scattering angle Θ. Equation 2 can be solved as

(3)

Thus, for the cross section dσ/d differential in the energy one derives

(4)

which is plotted in Figure 6. The figure shows that the largest part of the photon energy during Compton scattering is transferred onto electrons. Thus, the energy of photons is not shifted down slowly as in the energy range below 1 MeV [Celestin and Pasko, 2012]. As the number of new electrons produced through Compton scattering is small compared to the number of ambient electrons, we do not trace them in the rest of the section.

3.1.2 Generation of Hadrons

Photons also produce neutrons and protons in photonuclear reactions [Fuller, 1985] where a photon γ is absorbed by the nucleus of a molecule and a neutron n or a proton p is emitted:

(5)

(6)

Here Z is the atomic number, and M is the rest mass of atom A; we note here that the emission of protons produces atoms B with atomic number Z − 1 changing the composition of the ambient gas. In our model, we only use molecules of and as targets as the percentage of other nitrogen or oxygen isotopes in air is negligible. Since the binding energy E_bind of a nucleon is approximately 7.4 MeV for nitrogen and 8.0 MeV for oxygen (calculated with the Bethe Weizsäcker equation [Weizsäcker, 1935]), photons need tens of MeV to produce hadrons. The kinetic energy E_kin of the emitted hadron is

(7)

if the nucleus is left behind in its ground state. Hadrons are emitted isotropically; we neglect the motion of the residual nuclei as their rest mass is much higher than the rest mass of a neutron or a proton.

Figure 7a shows the total cross section for the photoproduction of hadrons from air molecules. It shows that the photoproduction of hadrons is most efficient for photon energies between 20 MeV and 25 MeV. Figure 7b shows the ratio of the number of produced hadrons to the number of produced positrons; it shows that in this energy range, the production of hadrons is 1 to 2 orders of magnitude less than the production of positrons.

3.1.3 Production and Motion of Positrons

We sample the total positron energy E₊ and positron direction Θ₊ relative to the direction of the incident photon, using the differential cross section of Bethe and Heitler in the integrated form of Köhn and Ebert [2014a]. We use the same elastic scattering, ionization, and Bremsstrahlung cross sections as for electrons. This is feasible since cross sections for electrons and positrons are similar for kinetic energies above 1 MeV [Agostinelli et al., 2003; Kothari and Joshipura, 2011]. We have also included the annihilation of positrons at shell electrons using analytic equations of Greiner and Reinhardt [1995]. Figure 8 shows the probability for annihilation of positrons at shell electrons of air molecules as a function of altitude for different positron energies when the positrons move from 16 km altitude upward excluding all other scattering processes. It shows that the probability is smaller than 15% for positron energies of 1 MeV and decreases rapidly with increasing positron energy.

3.1.4 Influence of the Geomagnetic Field

In order to estimate the influence of the geomagnetic field on relativistic electrons or positrons, we have to compare the gyration frequency in a magnetic field B with the collision frequency of electrons. For energies below 1 keV, these frequencies were compared in Ebert et al. [2010]. In general, the gyration frequency ν_G is

(8)

where we use the relativistic expression m(E_kin) = (E_kin + m₀c²)/c² for an electron with kinetic energy E_kin and e₀ and m₀ are the charge and the rest mass of an electron. The geomagnetic field is approximately 3 · 10⁻⁵ T at the equator up to altitudes of approximately 300 km. Note that ν_G is not constant but decreases with increasing electron energy because of the energy dependence of the electron mass m. The classical approximation ν_clas = e₀B/(2π·m₀) where the electron mass is taken as constant is also plotted in Figure 9. It shows that the gyration frequency 8 starts to deviate from the classical value for energies above 10 keV; for 1 MeV the gyration frequency is only one quarter of the classical value. The collision frequency ν_C is

(9)

where σ_tot(E_kin) is the total cross section as a function of E_kin, and where n_B(z) is the gas density as a function of altitude.

Figure 9 shows the comparison of the gyration frequency 8 and the collision frequency 9 for different altitudes for electron energies above 1 keV. It shows that for 16 km or 50 km and for energies between 1 keV and 100 MeV, the collision frequency is higher than the gyration frequency. Thus, the geomagnetic field is negligible. For approximately 120 km electrons with energies of ≈1 MeV start to feel the influence of the geomagnetic field. For approximately 140 km altitude the collision frequency is smaller than the gyration frequency for energies below 40 MeV. Hence, for altitudes between 120 km and 140 km, electrons and positrons with energies between 1 MeV and 40 MeV start to gyrate around the geomagnetic field lines. In our simulations we do not take the geomagnetic field into account. For altitudes below 120 km altitude and electron or positron energies above 1 MeV, we have just shown that the geomagnetic field is negligible; for altitudes substantially above 120 km beams of electrons and positrons follow the geomagnetic field lines, and their spatial distribution in the planes orthogonal to the field lines does not change any more.

3.3.1 Initial Condition

The solid line in Figure 4 shows that the number n_γ of photons with energy E_γ can be fitted well with

(11)

for photon energies above 1 MeV.

As shown in Figure 2, pair production and hadron production become relevant for energies above approximately 10 MeV. Since measurements [Briggs et al., 2010; Marisaldi et al., 2010; Tavani et al., 2011] have shown that TGFs can have energies of up to 40 MeV, we use the distribution 11 in the energy range from 5 MeV to 40 MeV as an initial condition and populate it with approximately 1.2 million photons. As the photons in our simulation are produced within 24 ns, thus within some meters and without much spatial separation, we initiate the photon beam at one single point at 16 km altitude at time zero. We use a monodirectional beam because most high-energy electrons move in forward direction and because the photons with energies above several MeV are emitted in forward direction. We trace the photon beam and its particle production for 1 ms which corresponds to a distance of approximately 300 km, and we use the barometric formula for the air density as a function of altitude.

3.3.2 The Energy and the Temporal Evolution of Positrons and the Energy of Hadrons

The photon beam 11 propagates upward with the scattering processes described earlier, and we have calculated the position and energy of photons and positrons after 0.5 ms. Figure 11 shows the position of photons (black) with energies above 5 MeV and of positrons (red) after 0.5 ms. As explained in 21, Figure 10 shows that a collision of a photon with an air molecule is very likely between 16 km and 20 km altitude. Thus, almost all photons have either disappeared due to the production of positrons or hadrons or lost so much energy through Compton scattering that their energy is less than 5 MeV and that they are removed from the pool of simulated photons. In contrast, Celestin and Pasko [2012] investigated photons with energies below 1 MeV where photons lose a smaller fraction of their energy through Compton scattering. Furthermore, Figure 11 shows that the positron beam moves with almost the speed of light.

Figure 12 shows the position and energies of all positrons including the positron scattering processes elastic scattering, ionization and Bremsstrahlung production. Figures 12a and 12b show the positron beam after 50 μs and 0.5 ms where the color denotes their kinetic energy. Figure 12c explicitly shows that positrons with energies above 20 MeV are in the front part of the positron beam while positrons with energies below 5 MeV are located rather at the end. Figure 12d shows the energy distribution of positrons after 1 μs, 50 μs, and 0.5 ms. It has a clear maximum at approximately 5 MeV. The shape of the distribution does not change considerably in time but only the total positron number.

Figure 13 shows the energy distributions of neutrons (Figure 13a) and protons (Figure 13b) after 10 μs, i.e., briefly after they have been produced, since we do not trace them through air. In both cases there are distinct maxima and minima due to the discrete structure of the photonuclear cross sections as shown in Figure 7. For neutrons the energies range from 4 MeV to 24 MeV; protons even have energies up to 33 MeV.

Our calculations are consistent with results of Babich [2007] and extend it. Babich [2007] calculates photon and neutron fluxes from an upward propagating atmospheric discharge with the help of cross sections and rate coefficients. He determines the mean energies of neutrons to be approximately 10 MeV which is consistent with the energy distribution of neutrons in Figure 13. Proton generation, however, has not been predicted before.