Production mechanisms of leptons, photons, and hadrons and their possible feedback close to lightning leaders
Abstract
It has been discussed that lightning flashes emit high-energy electrons, positrons, photons, and neutrons with single energies of several tens of MeV. In the first part of this paper we study the absorption of neutron beams in the atmosphere. We initiate neutron beams of initial energies of 350 keV, 10 MeV, and 20 MeV at source altitudes of 4 km and 16 km upward and downward and see that in all these cases neutrons reach ground altitudes and that the cross-section areas extend to several km2. We estimate that for terrestrial gamma-ray flashes approximately between 10 and 2000 neutrons per ms and m2 are possibly detectable at ground, at 6 km, or at 500 km altitude. In the second part of the paper we discuss a feedback model involving the generation and motion of electrons, positrons, neutrons, protons, and photons close to the vicinity of lightning leaders. In contrast to other feedback models, we do not consider large-scale thundercloud fields but enhanced fields of lightning leaders. We launch different photon and electron beams upward at 4 km altitude. We present the spatial and energy distribution of leptons, hadrons, and photons after different times and see that leptons, hadrons, and photons with energies of at least 40 MeV are produced. Because of their high rest mass hadrons are measurable on a longer time scale than leptons and photons. The feedback mechanism together with the field enhancement by lightning leaders yields particle energies even above 40 MeV measurable at satellite altitudes.
Key Points
- A terrestrial gamma-ray flash between 10 and 2000 neutrons per ms and m2 is possibly detectable at 0 km, 6 km, or 500 km altitude
- Combining the feedback mechanism with the electric fields of lightning leaders can explain beams of leptons and photons above 40 MeV
- The cross-section area of neutron beams between 350 keV and 20 MeV initiated at 4 km and 16 km altitude extends to several km2
1 Introduction
The last 30 years were marked with the discovery of several atmospheric phenomena, such as terrestrial gamma-ray flashes (TGFs) first observed by Fishman et al. [1994] as well as beams of electrons and positrons [Dwyer et al., 2008; Briggs et al., 2011]. There are various laboratory experiments confirming the production of high-energy electrons and photons correlated to electric discharges [Babich et al., 1990; Babich, 2003, 2005; Nguyen et al., 2008; Rahman et al., 2008; March and Montanyá, 2010; Shao et al., 2011; Kochkin et al., 2012, 2014] as well as satellite measurements supporting the idea of producing MeV particles correlated to lightning events [Smith et al., 2005; Briggs et al., 2010; Marisaldi et al., 2010; Tavani et al., 2011].
In particular Shah et al. [1985] performed the first measurements of neutrons related with thunderstorm and lightning activity. These measurements were done in the regime of 2.45 MeV which is the upper energy limit using detectors based on BF3. They estimated a production of 107–1010 neutrons per lightning discharge, but they did not make a hypothesis on the production mechanism of these particles. Further measurements were performed [Shyam and Kaushik, 1999; Bratolyubova-Tsulukidze, 2004; Tsuchiya et al., 2012; Kozlov et al., 2013, and references therein] enhancing the present knowledge about these phenomena. However, Babich et al. [2013a, 2013b, 2014] have shown that the counters in these experiments have been sensitive also to other ionizing emissions. More reliable observations were performed by Chilingarian et al. [2010, 2012].
Several simulations were performed in order to study the characteristics of these observations. Babich [2007] performed a theoretical analysis of enhancements in the neutron flux observed in association with lightning discharges. He showed that nuclear fusion in the lightning channel or related to thunderstorm electrical activity in general is not a probable mechanism for the neutron fluxes measured by Shah et al. [1985], Shyam and Kaushik [1999], Kuzhewskij [2004], and other authors. Instead, they propose that relativistic runaway electrons play an important role in atmospheric breakdown and that neutron bursts may be created by photonuclear reactions associated with this process. They estimate the number of relativistic electrons generated by upper atmospheric discharges, such as sprites and other transient luminous events, to be 1017 and that downward propagating neutrons generated by these discharges would get strongly attenuated, such that they would not be measurable by ground detectors. They also estimated that lightning leader process would yield 1013–1015 neutrons produced via photonuclear reactions and further suggested that this production is in better agreement with the experimental observations.
Grigoriev et al. [2010] developed simulations of neutron propagation through the atmosphere up to orbital altitudes in order to compare with the observations of Bratolyubova-Tsulukidze [2004]. They used GEANT4 to simulate the neutron propagation in order to estimate the possibility to obtain the neutron flux at orbital altitudes generated by photonuclear reactions. However, they did not provide any further analysis for beams of different energy or source altitude. Their results suggested that it is unlikely that these spaceborne observations are connected with TGF phenomena.
Carlson et al. [2010a], also using the GEANT4 software package, presented results of neutrons exiting a point-like source in three different distributions. They simulated the neutrons exiting isotropically from a point of 10 km initial altitude, an upward neutron beam with an initial altitude of 20 km and distributed over a 20° cone aperture, and a downward neutron beam with an initial altitude of 5 km, also distributed over a 20° cone aperture. Their simulations agree with ground observations such as the observations performed by Shah et al. but were not consistent with the spaceborne observations of Bratolyubova-Tsulukidze [2004].
To remove neutrons from air molecules by photonuclear reactions, photons with energies higher than 10 MeV are needed. The existence of such photons is correlated with the occurrence of terrestrial gamma-ray flashes (TGFs) and terrestrial electron beams (TEBs) [Dwyer et al., 2008]. There are two models explaining electrons and photons in the MeV range within the atmosphere: the Relativistic Run-Away Electron Avalanche (RREA) model [Wilson, 1925; Gurevich et al., 1992; Dwyer, 2003, 2012; Babich et al., 2012; Gurevich and Karashtin, 2013] and the acceleration of electrons in the vicinity of lightning leaders [Moss et al., 2006; Carlson et al., 2010b; Babich et al., 2011; Celestin and Pasko, 2011; Babich, 2013, 2015; Xu et al., 2012a; Köhn and Ebert, 2015].
Both models predict energies of photons and electrons of up to several tens of MeV. Besides the existence of high-energy electrons and photons, also positron beams have been measured [Briggs et al., 2011].
Dwyer [2003] was the first one to suggest a feedback model which includes the production of electron and positron pairs as well as the production of Bremsstrahlung photons by these leptons (in this article: electrons and positrons). Tracing all these particles in rather uniform electric fields, the feedback is capable of enhancing electric fields more than conventional models do.
In previous work [Köhn et al., 2014; Köhn and Ebert, 2015], we adopted the model by Xu et al. [2012a] and modeled the production of photons in the vicinity of a negative, stationary lightning leader. Subsequently, we modeled the production of beams of positrons, neutrons, and protons associated with TGFs using a Monte Carlo approach [Köhn and Ebert, 2015] in a two-step manner: first calculating the energy distribution of photons close to the lightning leader and second calculating the generation of positrons and hadrons (here neutrons and protons) from these photons. With this model approach, we have observed that all these particle species can reach energies of up to tens of MeV.
Although there is already a variety of work dealing with the feedback mechanism [Dwyer, 2007; Kutsyk et al., 2011; Skeltved et al., 2014], there is still a lack of feedback models which include all relevant particle species and interactions simultaneously rather than simulating first electrons then other species like photons, positrons, or hadrons. This is of importance since photons can create electrons and positrons which vice versa can produce Bremsstrahlung photons. Additionally, positrons can produce electrons via impact ionization. Protons are produced not only by Bremsstrahlung photons but also by neutrons which themselves are produced by photons. By tracing all species simultaneously including their collisions with air molecules, we include the production of particles even at higher altitudes; this effect is neglected if the motion of particles is modeled in a stepwise manner.
A full model should combine the generation and motion of leptons, photons, and hadrons in various field configurations and additionally include the above mentioned dependencies among different species. In the scope of this paper we introduce a full three-dimensional feedback model.
1.1 Organization of the Paper
This paper is organized in two parts.
In section 2 we introduce our three-dimensional, relativistic Monte Carlo feedback model and put emphasis on the interaction of neutrons with air molecules. We will discuss the main neutron dissipation processes. Moreover, we will introduce the concept of atmospheric depths which gives a mean to estimate the absorption and energy loss of particle beams traveling at the same air package but through different altitudes.
Finally, in section 3.1 we present simulation results of neutron beams with different initial energies originating at different altitudes and investigate the energy dissipation and how far they reach in the atmosphere.
In sections 3.2 and 3.3 we consider more general cases and analyze photon and electron beams of different energies. We explore the subsequent production and motion of electrons, positrons, protons, and neutrons and take into account the electromagnetic feedback of photons, positrons, and electrons. For these cases we then present the spatial and energy distributions after different time steps.
Our conclusions are presented in section 4.
2 Model
2.1 Cross Sections and the Monte Carlo Particle Code for Leptons, Photons, and Hadrons
We model the motion and production of leptons (electrons and positrons), neutrons, and photons as well as the production of protons in air, i.e., in a mixture of N2 (78.12%), O2 (20.95%), and Ar (0.93%), with a three-dimensional, relativistic Monte Carlo code where we explicitly model the collision of electrons, positrons, photons, and neutrons with air molecules.
For electrons we include elastic scattering off air nuclei [LXCat Database, 2014; Jacob, 1973; Phelps and Pitchford, 1985; Phelps, 1985; Fiala et al., 1994; Vahedi and Surendra, 1995], electron impact ionization [Yong-Kim and Santos, 2000], attachment to molecular oxygen [Luque and Gordillo-Vázquez, 2011; Pancheshnyi et al., 2008; Lawton and Phelps, 1978], molecular excitations [LXCat Database, 2014; Lawton and Phelps, 1978], electron-electron Bremsstrahlung [Köhn et al., 2014], and electron-nucleus Bremsstrahlung [Bethe and Heitler, 1934; Köhn and Ebert, 2014a]. For positrons we include the same processes plus the annihilation of positrons at shell electrons of air molecules. This is justified since positrons and electrons behave alike [Agostinelli et al., 2003Kothari and Joshipura, 2011].
For photons we include photoionization [EPDL Database, 1997], Compton scattering [Greiner and Reinhardt, 1995; Peskin and Schroeder, 1995] off shell electrons, the production of electron-positron pairs [Köhn and Ebert, 2014a], and the production of neutrons and protons [Fuller, 1985]. Through the removal of neutrons from 
N or 
O by photons, the radioisotopes 
N or 
O are produced. However, as we have shown in [Köhn and Ebert, 2015], this process is not significant enough to produce radioisotopes in large amounts.
For neutrons we include elastic scattering, excitations of air nuclei, and charge exchange reactions, i.e., A(n, p)B, where A is the nucleus of an air atom (target nucleus) and B is the nucleus formed by exchanging a proton by a neutron in the target nucleus. There we use different cross sections for nitrogen, oxygen and argon [Mendoza et al., 2012, 2014]. The charge exchange process can eventually produce radioisotopes, e.g., as through 
N(n, p) 
C. However, this process is a third-order process (after the production of photons and the subsequent production of neutrons), and thus, it is negligible.
By implementing all the reactions mentioned above, we have automatically taken into account the feedback of electrons and positrons producing new photons through the Bremsstrahlung process and through positron annihilation at shell electrons, as well as of photons producing leptons through Compton scattering, photoionization, and pair production. The treatment of hadrons is not relevant for the feedback itself, but we have included them to complete the microphysical part as much as possible. However, once neutrons are liberated from air nuclei they can produce protons through charge exchange reactions and, as such, enhance the proton signal.
We note here that this is the same microphysics as in [Dwyer, 2003] where new RREAs are initiated by electrons produced by photons or positrons. However, in our model we do not look into RREAs in homogeneous fields. Hence, our feedback terminology differs slightly from the one used in Dwyer [2003] since we do not look into the feedback of photons and positrons on RREAs but rather into the feedback of one species on another species.
The Monte Carlo code consists of two alternating steps: Between two steps the position and the velocity of the particles are updated according to the relativistic equations of motion in a given electromagnetic field. Afterward, we draw random numbers to determine if a collision, and if so, which collision takes place. For all species we use the same Monte Carlo algorithm as used by Li et al. [2012] fed with the above mentioned cross sections for different species. The probability of a particle colliding with an air molecule is proportional to the air density which we control with the barometric formula. Details of the implementation of the electromagnetic processes can be found in [Köhn and Ebert, 2014b, 2015].
Parts of this Monte Carlo code were used to calculate the energy distributions of photons and positrons in the vicinity of a negative stepped lightning leader [Köhn and Ebert, 2015] and the production of X-rays from meter-long discharges [Kochkin et al., 2016]. We already validated our code up to energies of 1 MeV in [Kochkin et al., 2016] by comparing the National Institute of Standards and Technology (NIST) data [Berger et al., 1998] with our simulation results for the average distance ν for full energy loss of an electron with initial energy E0. The NIST data are the so-called CSDA (continuous-slowing-down approximation) range which is “a very close approximation to the average path length traveled by a charged particle as it slows down to rest, calculated in the continuous slowing-down approximation. In this approximation, the rate of energy loss at every point along the track is assumed to be equal to the total stopping power. Energy loss fluctuations are neglected. The CSDA range is obtained by integrating the reciprocal of the total stopping power with respect to energy.” (NIST CSDA, http://physics.nist.gov/PhysRefData/Star/Text/appendix.html) In the scope of this paper we extended this comparison up to 50 MeV, shown in Figure 1a. Additionally, Figure 1c shows the relative error (νNIST−νSimulation)/νSimulation indicating good agreement between the data and the simulations with a relative error of 
0.2 except for a small deviation for very small and very large energies. However, for intermediate energies where most of the electrons are in the energy distribution, the agreement is much better. As mentioned above, for positrons we use the same cross sections as for electrons in addition to the annihilation of positrons at shell electrons.
Figure 1b compares the mean-free path Λ of photons in air tabulated in the NIST database [Berger et al., 1998] and calculated with our Monte Carlo code; Figure 1d) compares the relative difference (ΛNIST−ΛSimulation)/ΛSimulation. There are the deviations between the NIST data and our simulation results 
0.4 for photon energies below 1 MeV; for photon energies above 1 MeV, however, the relative error tends to zero. This is particularly the energy range we are interested in for the production of electron-positron pairs and of hadrons.
2.2 Neutron Cross Sections
In previous work, e.g., Dwyer [2007] and Köhn and Ebert [2015], cross sections for electrons, positrons, and photons have carefully been investigated and documented. Here we revise cross sections for neutrons moving through the atmosphere.
(1)
km) where n0≈2.6·1025 m−3 is the air density at standard temperature and pressure. We have made a similar analysis for photons in [Köhn and Ebert, 2015]. If 1/Nint(r1,r2) < σ for a process with cross section σ , a collision of a neutron with an air molecule is very likely. Figure 2b shows that for neutron energies above 2 MeV inelastic scattering is very probable for the motion of neutrons between 16 km and 500 km or between 16 km and 0 km; neutrons traveling between 16 km and 19 km scatter rather elastically than inelastically.
, [Mendoza et al., 2012, 2014] where Θ is the scattering angle. Once Θ is known, we can use the conservation of energy and momentum to determine the kinetic energy Ekin,n,f of the neutron in the final state; it is
(2)
(3)
(4)
(5)
For energies above 10 MeV, the cross section of inelastic scattering approaches the cross section of elastic scattering. Inelastic scattering of neutrons can either be the excitation of air molecules, their nuclei, or the production of charged nucleons through the (n, p) charge exchange reaction. Figure 4 shows the average energy loss ΔE (solid line) for excitations as a function of the incident neutron energy. It shows that the energy loss grows linearly with the energy of the incident neutron. Additionally, the dotted line shows the fraction of the energy loss to the incident neutron energy Ekin,n,i. On average more than 60% of the energy, much more than through elastic scattering, is lost through excitations.
The propagation of neutrons through the atmosphere is an interplay of inelastic and elastic scatterings. After neutrons have lost energy through inelastic scattering, the cross section for elastic scattering increases and the average scattering angle becomes larger. Thus, the number of elastic scatterings per time unit increases, and subsequently, the energy loss through elastic scattering becomes more and more significant.
Section 3 not only presents the attenuation of neutrons in air but also will show that the implementation of the above mentioned cross sections and algorithms yields attenuation lengths comparable to those mentioned in literature. For further validation we compare with data provided by Nakamura and Kosako [1981]. They calculate the total neutron flux r·Φ(r) of an almost monoenergetic upward directed neutron beam at distance r from the source for different initial energy groups where the flux Φ = Nn/(4πr2) is defined as the number of neutrons Nn passing through a sphere with radius r. Figure 5 shows the comparison of the flux of a neutron beam for different distances calculated with our simulation tool with the results presented in Figure 8 of [Nakamura and Kosako, 1981]. In [Nakamura and Kosako, 1981] they do not use distinct initial energies but rather energy groups. We compare a 1.5 MeV beam with a beam in energy group 14 (1–2.02 MeV) and a 150 keV beam with energy group 16 (67.4–247 keV). We observe the same exponential decrease and a good agreement between our results and the ones by Nakamura and Kosako [1981]. We see some deviations between our simulations and those of Nakamura and Kosako [1981] which might be related to the fact that they use energy groups rather than distinct initial energies.
2.3 Altitude and Atmospheric Depth
(6)
3 Results
We have simulated the motion of different initial monoenergetic particle beams of different species with different initial energies E0 from different altitudes H0, i.e., air package ΔX = 0, in the atmosphere and calculated the temporal evolution of the spatial and energy distributions of all relevant particles. We have chosen the sign of ΔX such that ΔX < 0 refers to altitudes in the initial direction of the particle beam.
The following three subsections establish the feedback model in a piecewise manner. In section 3.1 we look into neutrons, in section 3.2 we look into photons which eventually produce neutrons, and in section 3.3 we look into beams of electrons which can produce photons subsequently liberating neutrons from air nuclei.
The duration of the simulations of neutron beams was approximately 2–3 weeks on an “Intel(R) Xeon(R) CPU E3-1270 v3 at 3.50 GHz;” the simulations of the particle beams in the vicinity of lightning leaders took approximately 6 months on an “Intel(R) Xeon(R) CPU E5-2650L v2 at 2.60 GHz.”
3.1 The Motion of Monoenergetic Neutron Beams Through the Atmosphere
We simulated the motion of monoenergetic neutron beams originating at different altitudes upward and downward; all beams are initiated with 600 neutrons. We note here that we only regard some particular cases to cover the investigation of neutron beams in the keV and MeV range which is the energy range of neutrons as being produced by TGFs [Köhn and Ebert, 2015]. Hence, we get a broad insight into the properties of neutron beams in the atmosphere.
Figure 6 shows the spatial distribution of neutrons after 1 ms for initial energies of 350 keV (a), 10 MeV (b–d), and 20 MeV (e) for beams from 4 (a, b, and e) and 16 km (c and d) altitude, upward (a and d) and downward (b, c, and e). Figure 7 shows the energy distribution after 0.25 ms, 0.5 ms, 0.75 ms, and 1 ms for the same cases. Neutrons with energies of 350 keV are rather scattered isotropically and travel less than 10 km in air; in the upward directed beam (Figure 7a) some of the neutrons are backscattered and reach the ground.
The color of the points as well as Figure 7a show that there is a maximum in the energy distribution at 100 keV and also that there are still neutrons with 350 keV, thus keeping all their initial energy. As Figure 3 shows, inelastic scattering is not dominant for energies below 500 keV, and the only energy loss is through elastic scattering. Due to the randomness of neutrons scattering at air molecules, some neutrons will collide more frequently and lose more energy whereas a few of them will not collide too often and thus keep their energy.
Figures 6b, 6e, 7b, and 7e show that beams of 10 MeV and 20 MeV emitted from 4 km downward behave similarly. After 1 ms neutrons populate altitudes from 0 km up to 8 km; because of the high probability of inelastic scattering for energies of 10 MeV and 20 MeV the neutrons lose a significant ratio of their initial energy, and hence, after 1 ms most neutrons have energies below 80 keV. The situation changes if the neutron beam was emitted at 16 km altitude (Figures 6c and 6d). For both upward and downward emitted beams there are only a few neutrons reaching the ground. There is a cloud of neutrons close to the source altitude and also neutrons moving upward. Figures 7c and 7d show that for the downward emitted beam most neutrons have energies below 1 MeV whereas for the upward emitted beam there is a significant number of neutrons with energies of approximately 10 MeV. While moving upward, these neutrons have not undergone a considerable number of collisions in the early stage of the beam evolution, and once the neutrons have passed a certain altitude, the air density has thinned out; thus, collisions become less probable. For the downward emitted beam neutrons, they have first to backscatter before they can move upward; since neutrons lose most of their energy when they are backscattered, the neutron energies are smaller than the ones for the upward emitted beam.
Figure 8 shows the neutron number N normalized to the initial neutron number N0 passing through altitude z. Figure 8a shows an overview of all simulated cases. It shows that the neutron number decreases as the neutrons travel through the atmosphere. This is due to inelastic collisions, and for all simulations there is a small probability of detecting neutrons on the ground according to Shah et al. [1985], Shyam and Kaushik [1999], Bratolyubova-Tsulukidze [2004], Tsuchiya and other [2012], and Kozlov et al. [2013]. It also shows that for the same initial energy and altitude showers behave alike no matter whether they move upward or downward as long as they pass through the same air package ΔX. Note here that equation 6 is only valid for altitudes below 112.8 km; hence, the green line representing an upward emitted beam of 10 MeV neutrons starting at an altitude of 16 km is interrupted abruptly. ΔX = 0 means no air layer at all and denotes the initial altitude of the neutron beam. Since, we measure the flux of neutrons relative to ΔX = 0, N/N0 drops down to 0 at ΔX = 0.
The attenuation length λatt (g/cm2) is defined as the atmospheric depth where the neutron number has dropped to 1/e of the initial neutron number (in the initial direction). For 10 MeV neutrons, it is λatt≈74 g/cm2; for 20 MeV neutrons it is λatt≈84 g/cm2 which is in agreement with the values presented by Shibata [1994]; Koi et al. [1999].
Figures 8b–8f show details for all simulated cases. In each panel for different initial energies, different lines show the neutron number as a function of altitude/atmospheric depth for above different energies Ekin,n; in these panels Ekin,n≥0 means that all neutrons are taken into account, hence, showing the same line as in Figure 8a. Figures 8b–8f explicitly show that most neutrons move in their initial direction. Comparing Figure 8a with the others shows that for a 350 keV beam more neutrons move in the opposite direction than for beams of 10 MeV or above. Figures 8d and 8e show that for a 10 MeV beam emitted at 16 km altitude upward or downward, hardly any neutrons reach the ground but move up to several tens of km. The different lines in Figures 8b–8f show that, regardless of the initial energy, there is a vanishing flux of high energy neutrons and that neutrons reach the ground after having lost a significant fraction of their initial energy. Because of the small neutron scattering cross section for high neutron energies, neutrons with higher energies travel further before losing their energy than low-energy neutrons which leads to the staircase structure in Figures 8b–8f. The staircase structure is due to discretized output of our Monte Carlo simulation every 10 μs. As such, especially for high neutron energies, we cannot plot the decrease of the neutron number accurately enough and sudden false stairs occur.
Figure 9 shows the normalized cross-section area 
as a function of H or atmospheric layer ΔX. N(H) is the total neutron number at altitude H and xi(H) as well as yi(H) are the coordinates of the ith neutron at altitude H parallel to the xy plane. Figure 8a shows an overview over all simulated cases whereas Figures 8b–8f show details for each simulation. Figure 8a shows that for all beams from 4 km downward and upward the neutron beam covers an area of approximately 1 km2 at H≈0. At the ground the 10 MeV beam from 16 km downward covers an area of approximately 10 km2 whereas the 10 MeV beam from 16 km upward covers an area of approximately 50 km2 because on their way downward these neutrons scatter more than those initiated at 4 km altitude.
as a function of atmospheric layer ΔX and equivalently of the altitude H. xi and yi are the positions of the ith neutron (out of N) parallel to the xy plane. The solid line at ΔX = 0 indicates the source altitude of the neutron beam. (a) All simulated cases; (b–f) 〈A〉 for different threshold energies Ekin,n for different initial altitudes and different initial energies.
Figures 6 and 9 also show that for all simulated cases the neutron beams are much more widened for high altitudes in contrast to small altitudes. This is very well illustrated in Figures 6d and 9e. Close to the source location of the beam, neutrons scatter at air molecules, and hence, their direction changes; once the neutrons reach a region with a much smaller air density, they do not scatter anymore and continue moving in the direction they obtained at lower altitudes. Neutrons moving downward are confined because they keep on scattering at air molecules.
For all considered cases, Table 1 displays the fluence Ψ = Nn/(Δt·A) of neutrons at different altitudes, where Nn is the number of neutrons passing through the cross-section area A within the time interval Δt; a small fluence is equivalent to a small neutron number or a large cross-section area. Table 1 shows the fluence at ground (0 km), at typical balloon experiment altitudes (6 km), and at typical satellite altitudes (500 km); for those cases where our simulations stopped before neutrons reach the considered altitude or where all of the test neutrons are absorbed beforehand, we present upper bounds of the fluence. For almost all cases, the fluence lies between 10−11 and 10−9 per m2 and ms. The only outliers appear for the 10 MeV neutron beam initiated at 16 km altitude downward; in that case the flux at ground and at satellite altitude is below approximately 10−13 m−2 ms−1.
| Initial Condition | H (km) | Nn/(Δt·A) (1/(m2 ms)) | Nn,TGF/Δt (1/ms) |
|---|---|---|---|
| 350 keV/4 km/up | 0 | <10−9 | <1667 |
| 350 keV/4 km/up | 6 | 5·10−9 | 8333 |
| 350 keV/4 km/up | 500 | <10−9 | <1667 |
| 10 MeV/4 km/down | 0 | 5·10−9 | 8333 |
| 10 MeV/4 km/down | 6 | 5·10−10 | 833 |
| 10 MeV/4 km/down | 500 | <10−10 | <167 |
| 10 MeV/16 km/down | 0 | <10−13 | <0.17 |
| 10 MeV/16 km/down | 6 | 3·10−10 | 500 |
| 10 MeV/16 km/down | 500 | <2·10−13 | <0.33 |
| 10 MeV/16 km/up | 0 | <10−11 | <17 |
| 10 MeV/16 km/up | 6 | 10−10 | 167 |
| 10 MeV/16 km/up | 500 | <2·10−11 | <33 |
| 20 MeV/4 km/down | 0 | 10−9 | 1667 |
| 20 MeV/4 km/down | 6 | 10−10 | 167 |
| 20 MeV/4 km/down | 500 | <10−11 | <17 |
The total number of neutrons detectable at altitude H per time interval related to a TGF event can be estimated through Ψ·Nn,TGF/600·Adet where Nn,TGF is the number of neutrons produced during one TGF flash, 600 the initial number of neutrons in our simulations and Adet the size of the detector. Here we choose Nn,TGF=1015 [Babich, 2007]; for the detector size we choose approximately 1 m2 for ground detectors [Chilingarian et al., 2010], for balloon experiments [Marshall T. C., et al., 1995], and for detectors in satellites [Meegan et al., 2009]. Table 1 shows that for most relevant cases the number of neutrons passing through a detector per millisecond lies between 10 and 2000 with some outliers for a 350 keV beam initiated upward (8333/ms) and for a 10 MeV beam initiated downward at 16 km (0.17/ms).
3.2 The Motion of Photons Through the Atmosphere and the Subsequent Production and Motion of Leptons and Hadrons
This section is devoted to investigate the dissipation of already existing photons in the vicinity of an ideally conducting lightning leader. The full ab initio calculation initiated with low-energy electrons will follow in the next section.
(7)
We investigate three different cases: no ambient field at standard temperature and pressure (STP); a uniform, constant electric field of −100 kV/cm at STP; and the field of a negative stepped lightning leader [Köhn and Ebert, 2015] of 1 km length and 1 cm curvature radius whose tip is at 4 km altitude, in an ambient field of −0.5 kV/cm, as used by Xu et al. [2012a]; Köhn and Ebert [2015]. Figure 10 shows the electric field in the vicinity of the lightning leader as calculated by Köhn and Ebert [2015]. There is a large field enhancement close to the leader tip which is enough to accelerate charged leptons into the run-away regime in the vicinity of the leader tip.
are used, and the upper leader tip lies at the origin of the coordinate system. The black level lines indicate fixed values of the electric field strength from 5 to 1000 kV/cm as indicated.
While propagating through the atmosphere, photons produce neutrons and protons as well as positrons and electrons producing new Bremsstrahlung photons. For the lightning leader configuration Figure 11 shows the detailed position and energy for all different species separately. We note here that we have not implemented collision processes for protons yet and, as such, on their path through the atmosphere they are only affected by electric fields. In future work we will elaborate their correct motion, too.
There is a distinct beam of photons upward. While moving upward, they produce leptons and hadrons. In the presence of any electric field the motion of high-energy electrons is directed against the electric field lines whereas positrons and protons move along the electric field. Since hadrons are emitted isotropically, they move through the atmosphere in such an isotropic manner where the motion of protons is superposed by the electric field. Figure 11 demonstrates very well the higher mass of neutrons and protons relative to those of leptons and photons. Because of that distinction, hadrons move much more slowly through the atmosphere than leptons with the same energy. Depending on the electric field, there is the build-up of different layers of species. Dwyer [2003] already calculated this separation for electrons and positrons. Figure 11 shows that in the presence of an electric field there is an additional separation due to the different charges of electrons on the one side and positrons as well as protons on the other side. Since in Figure 11 the electric field is not uniform, negatively charged electrons are accelerated not only upward but also sideward. Likewise, positively charged protons and positrons are accelerated not only downward but also sideward. Hence, charged particles are not beamed along the z axis.
Figure 12 shows the energy distribution of photons, leptons and hadrons after 1 μs (a and b) and after 0.5 ms (c and d) for energies above 5 MeV in the case of no electric field (a and c) and of the electric leader (b and d). Figures 12a and 12b show that in both cases after 1 μs there are leptons and hadrons with energies of up to approximately 30 MeV. Note that the fraction of photons converted to leptons and hadrons is small [EPDL Database, 1997; Köhn and Ebert, 2014a; Fuller, 1985], and thus, there are many fluctuations in the lepton and hadron numbers. Especially, for all simulated cases, there is a hole in the energy distribution of neutrons at approximately 11 MeV which we have already observed in previous work [Köhn and Ebert, 2015]; this hole is an artifact of the resonances to be found in the cross sections for neutron production (see Figure 7a in [Köhn and Ebert, 2015]). Figures 12c and 12d show that most photons have dissipated after 0.5 ms and that the number of leptons and hadrons has grown. In both cases there are even some rare events where we find hadrons with energies of up to 40 MeV. In the presence of an electric leader field, some few electrons have gained energies of up to 50 MeV. These electron energies are consistent with energies in recently observed electron beams [Dwyer et al., 2008]; in addition, Figure 12 shows that a typical photon beam as the initial state of a TGF can produce neutrons of up to several tens of MeV.
3.3 The Motion of Electrons in the Electric Field of a Negative Stepped Lightning Leader
We also simulated the motion of a monoenergetic electron beam with initial energy of 0.1 eV upward from 4 km altitude in the field of a leader of 1 km length and of 1 cm curvature at the tip, thus in the same leader field as in section 3.2. As we have seen in previous work [Köhn and Ebert, 2015] and in Figure 10 the enhancement of the electric field through the leader is large enough to overcome the friction of air and accelerate electrons into the run-away regime. As in [Köhn and Ebert, 2015] and in [Xu et al., 2012a], we use a simplified model where we do not take space charges into account. After they have gained enough energy, they produce photons through Bremsstrahlung which subsequently produce other species. Figure 13 shows the energy distributions after 28 ns and 0.5 ms. It shows that after 28 ns there are already electrons and photons with energies of up to 9 MeV and 16 MeV, respectively. After 0.5 ms electrons have reached energies of up to 50 MeV, and there are even positrons with energies of up to 50 MeV. We note here that, as in section 3.2, the positron spectrum has a maximum at approximately 5–10 MeV as has been calculated in [Köhn and Ebert, 2015] and been measured by Fermi. Figure 13 also shows that there are single hadrons with energies of tens of MeV.
After 0.5 ms leptons, photons and hadrons have reached altitudes of up to 150 km and reached energies up to 50 MeV. Because of the thinning of the air density as a function of altitude, most of these leptons will continue moving upward, following the geomagnetic field, and can be measured by satellites. Energies of particles in electron and positron beams as well as in TGFs are measured to reach to at least 40 MeV [Briggs et al., 2010; Marisaldi et al., 2010; Tavani et al., 2011] as is also in Figure 13b. Since in our model leptons can gain energies of above 40 MeV and since there is a small but not vanishing probability that these high-energy electrons might produce Bremsstrahlung photons, even at higher altitudes, this sheds light on the fact that photons with energies above 40 MeV, possibly even of up to 100 MeV might reach satellite altitudes (M. Marisaldi, personal communication, 2015); this effect is not observable if we trace electrons only in the beginning, i.e., some meters in the vicinity of the leader tip although the leader field is still present several kilometers ahead of the leader tip, and not at the later stage of the simulation or if we do not model the production of electrons by Bremsstrahlung photons at higher altitudes. We observe that these energies can be reached by initiating an electron beam at 4 km altitude. Additionally, as we have seen in section 3.1, there is also a small probability that neutrons which are produced by upward moving photons can be detected not only by satellite altitudes but also by ground-based instruments.
4 Conclusions and Outlook
In the course of this paper we investigated the motion of different beams of particles through the atmosphere.
The investigation of different neutron beams has shown that there is a small number of neutrons with energies of up to 20 MeV initiated upward and downward reaching the ground and also that upward initiated neutron beams with energies above 10 MeV can be detected by satellites. However, we have seen that the cross-section areas of these beams are in the order of 0.1–100 km2 which makes them harder to detect, especially for large cross-section areas. For an average terrestrial gamma-ray flash, producing 1015 neutrons, approximately 10 to 2000 neutrons per millisecond, depending on their initial energy and initial altitude, are detectable at ground, at 6 km or 500 km altitude. By analyzing cross-section data, we have seen that most of the energy loss originates from inelastic channels whereas elastic collisions do not contribute significantly to the energy loss.
From the investigation of different initial photon and electron beams in different electric field configurations we conclude that in all relevant cases there is the production of leptons (electrons and positrons), photons, and hadrons (protons and neutrons) with energies of up to several tens of MeV. Especially, in the case of a lightning leader of 1 km length, we have observed that electron or photon beams already initiated at 4 km altitude can produce these energies. However, we have seen that in electric fields high-energy electrons and photons move much faster ahead than protons, neutrons, and positrons. Positrons and protons are dragged into the opposite direction of the electrons because of their positive charge; even in the absence of any electric field, protons and neutrons move much more slowly than leptons because of their heavier rest mass. Hence, we expect hadrons to be measured on a much longer time scale than leptons and photons.
In future work also cross sections for protons should be implemented into our Monte Carlo code to trace them through the atmosphere correctly. Additionally, the effect of space charges should be included in future work to study their effect on the low-energy part of the particle spectrum. Furthermore, our Monte Carlo code can be used to study how a particle signal at ground or at satellite altitude changes if parameters like the initial particle energy, the relevant altitude or the electric field as well as the geometry of the leader are altered. This could potentially involve electron beams close to lightning leaders or monoenergetic neutron beams.
Acknowledgments
This work was supported by the European Science Foundation research networking project “Thunderstorm Effects on the Atmosphere-Ionosphere System” (TEA-IS). We thank Ute Ebert, Olaf Scholten, and Fernanda São Sabbas for helpful discussions. C.K. acknowledges financial support by STW-project 10757, where Stichting Technische Wetenschappen (STW) is part of The Netherlands' Organization for Scientific Research NWO. The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement 609405 (COFUNDPostdocDTU). G.D. acknowledges financial support by Brazilian institutions CAPES and CNPq. The simulation code can be obtained from C.K. ([email protected]) and will soon be published on the web.
state of O2 by low-energy electrons
