Runaway electrons from a ‘beam-bulk’ model of streamer: application to TGFs

The generation of x- and gamma-rays in atmospheric discharges has been studied intensively since the discovery of terrestrial gamma-ray flashes (TGFs) by the Compton gamma-ray Observatory in 1991. Emissions are bremsstrahlung from high energy particles accelerated in large scale atmospheric electric fields associated with thunderstorms. Whereas observations now are many, both from lightning and the laboratory, the phases of the discharge where emissions are generated are still debated and several processes for electron acceleration have been put forward by theorists. This paper address the electron acceleration in streamer region of lightning. We present the first ‘beam-bulk’ model of self-consistent streamer dynamics and electron acceleration. The model combines a Monte Carlo Collision code that simulates the high-energy electrons ( $?> > ?> 100 eV) and a fluid code that simulates the bulk of the low-energy electrons and ions. For a negative streamer discharge, we show how electrons are accelerated in the large electric field in the tip of the streamer and travel ahead of the streamer where they ionize the gas. In comparison to the results obtained with a classical fluid model for a negative streamer, the beam-bulk model predicts a decrease of the magnitude of the peak electric field and an increase of the streamer velocity. Furthermore, we show that a significant number of runaway electrons is lost by diffusion outside of the streamer tip. The results presented here do not yet include extra amplification nor acceleration far away from the streamer to explain the electron energies seen in TGFs. Still, in the light of those results, we emphasize that the production of runaway electrons from streamers needs to be simulated including the self-consistent feedback of runaways on the streamer. Simulations with a beam-bulk model may not only help to understand the fundamental atmospheric processes behind TGFs, but also pave the way for the interpretation of remote sensing of the most energetic discharges in the Earthʼs atmosphere and thus help to address their environmental impact.


Introduction
During the last 25 years interest in atmospheric electricity has increased because of the discovery of electromagnetic and radiative processes occuring in the atmosphere above thunderstorms. The first was transient luminous events (TLEs) (e.g. [2][3][4][5][6]), first observed in 1989 [1].
The second was sub-ms duration bursts of photons from thunderstorms with energies reaching above 20 MeV, the terrestrial gamma-ray flashes (TGFs), observed first in 1991 from the compton gamma-ray Observatory when above thunderstorms [7]. TGFs are bremsstrahlung from highenergy electrons accelerated in thunderstorm electric fields (e.g. [8]). The acceleration mechanism, and its relation with lightning and TLEs, is still debated. We know that electric fields exist on several spatial and temporal scales. On the smallest scale is the field of the streamer tip, then follows the field around a leader as in lightning or a blue jet, then the space charge field in the thunderstorm clouds and finally, the impulsive fields generated by lightning discharges reaching from the clouds to the ionosphere. The challenge is to identify a plausible combination of those fields that will create pulses of high-energy electrons which, in turn, generate photon pulses that match TGF observations.
The theoretical models, put forward since the discovery of TGFs, are built on Wilsonʼs proposal that the electrons accelerated in thunderstorm fields may reach the runaway regime [9]. Indeed, the frictional force, induced due to collisions with the atmospheric constituents, decreases with increasing electron energy, thereby allowing electrons to reach relativistic energies [9]. The addition of Møller scattering to the theory led to the relativistic runaway electron avalanche (RREA) mechanism which gives an amplification in the number of runaway electrons of up to 10 5 , assuming large-scale fields in thunderstorms and above [10]. The addition of pair production, creating antimatter in the form of positrons, and x-rays and their interactions with the atmosphere, led to the relativistic feedback mechanism which gives an extra amplification of up to 10 13 [11].
The mechanisms proposed have in common that they require seed electrons with energies in the runaway regime. It was proposed early, that these electrons may come from cosmic ray interactions with the atmosphere. More recent studies investigate, with particle models, if cold electrons can be accelerated into the runaway regime in the small-scale fields of streamers feeding a leader, and if the number of those runaway electrons is high enough to explain observations [12][13][14][15][16].
The interest for TGFs is still growing, not only on the theoretical basis but also for their environmental impact. The high energy particles (photon, electron positron and neutron) associated with TGF modify the radiation environment at low altitudes. The radiation dose predicted has been shown to be of importance for aviation safety [17][18][19].
In this work, we propose to use an hybrid approach coupling fluid and particle models to study electron acceleration in streamers, the effect of energetic electrons on streamer propagation and the production levels of relativistic electrons.
Recently, hybrid codes have been used to simulate streamer discharges in coupling fluid and particle models in different spatial regions, with the fluid description in the streamer plasma channel, where the electron density is large, and the particle description ahead of the streamer front (e.g. [21][22][23]). In this work, we present an alternative approach proposed by Belenguer and Boeuf [20]: the beam-bulk model. In this model the electrons, populating the computational domain, are not separated in space, but are rather divided in the energy space into two distinct groups: (a) the low energy electron population described by a fluid model, and (b) the high energy electrons followed by particles. This paper is organized as follows: first, we present details of the beam-bulk model. Then we present a simulation that illustrates the impact of runaway production on the streamer propagation, and its dependence on the altitude (the neutral density). Finally, we discuss the production mechanism and the production rate of runaway electrons.

The beam-bulk model of streamers
The beam-bulk model is derived from the pure particle model of Chanrion and Neubert [13,14], which couples an electrostatic particle-in-cell code (PIC) with a Monte Carlo Collision code (MCC), by replacing the PIC-MCC model for the low-energy electrons with a classical drift-diffusion fluid model (e.g. [24][25][26][28][29][30][31]). The high-energy particles are still followed using the PIC-MCC code which is extended in energy with the relativistic binary-encounter-Bethe electron impact ionization model from Celestin and Pasko [32].
In this work as a first step, the streamer is simulated using a 1.5D model in a point-to-plane geometry (e.g. [33][34][35][36][37]) shown in figure 1. With a 1.5D model, the discharge is considered as a cylinder with a finite radius and uniform radial distribution of charges and the evolution of the charged particle densities is solved in one spatial dimension along the direction of streamer propagation. It is important to note that the 1.5D model captures the main characteristics of streamer propagation, but extension of the beam-bulk model to 2D and 3D will be the subject of future work. In this work, the driftdiffusion equations are solved in the fluid approximation and the trajectories of energetic electrons are integrated in the particle approximation. The electron and ion densities arising from the energetic particle updates are calculated using a standard weighting scheme [13]. With a 1.5D model, the electric field is found by slicing the plasma cylinder net charge density into discs. Using a 3D analytical formulation for the axial electric field from a disc, the space charge electric field in the electrode gap is calculated by integrating the individual contributions of the disks and their image charges in the electrodes [33]. All simulations in this work are presented without pre-ionization or photoionization in order to emphasize the role of high-energy electrons. The discharges are initiated with a neutral plasma seed composed of Figure 1. The studied geometry. The cylinders represent the streamer space charge densities used to calculate the electric field. The streamer propagates from the point electrode on the right, maintained at a potential V, to the grounded planar electrode on the left. electrons and ions at rest with a density that follows a gaussian of peak value − 10 cm 16 3 and standard deviation 0.5 mm. The boundary in electron energy between the drift-diffusion and the PIC-MCC models is chosen to 100 eV. This energy is the one of an electron having approximately the velocity of the streamer. The electron mobility, ionization rate, attachment rate and diffusion coefficient of the driftdiffusion model have been tabulated from pre-computed swarm calculations using the PIC-MCC code in constant electric fields and for electrons with energies below 100 eV. In this way we have secured consistency between the two models. Computer particles from the PIC-MCC model that loose energy to below 100 eV are injected in the fluid model by calculating their density over the grid mesh. The electrons from the fluid model injected into the particle model are calculated from a table of distribution functions pre-computed in a constant electric field. Knowing the electric field and the electron density on the mesh, this allows to determine the probability that electrons reach energies above 100 eV. The coupling of electrons across the energy boundary is repeated at every time step.
The code allows for long streamer simulations by keeping the streamer head centered in the middle of the simulation after it has reached this position. This is implemented by shifting the grid backwards for one grid mesh every time the maximum value of the electric field in the head passes the middle of the simulation box.

Evaluation of the beam-bulk model
In the following, we compare the beam-bulk model with a pure drift-diffusionmodel. In the pure drift-diffusionmodel we use for consistency the electron mobility, ionization rate, attachment rate and diffusion coefficient tabulated the same way as for the beam-bulk model, except that we do it for all electron energies. We also evaluate the beam-bulk model relative to the neutral gas density n gas . Since there is no photoionization included in the models, the scaling is quite simple [26]. Time, t, and distance, x, scale as n 1 gas , and the electron density, n e as n gas 2 . The results presented in the following for different altitudes are scaled to sea-level density, unless otherwise noted. We present results on a negative streamer propagation for an electrode tip of radius μ 3.5 m, a gap of 1 cm, and an applied voltage of −25 kV. The streamer radius is fixed to 0.5 mm at ground pressure (see [38] and ref therein). The simulation time step is − 10 ns 4 and the mesh size is μ 5 m. The axial profiles of electron density and electric field at = t 0.8 ns are shown in figure 2. Results of the pure drift-diffusion model is shown on the two top panels. When scaled, it is independent of the altitude (n gas ). Below are shown results of the beam-bulk model for three different altitudes: 0 km, 10 km and 20 km. The left panels show the full simulation domain and the right panels a close up of the streamer head. The electron density is shown in blue, the density of electrons above 100 eV is in red and the magnitude of the electric field, normalized by the local conventional breakdown field E k , is shown in green.
In figure 2, the streamer propagates from the right towards the left. We see that the high-energy electrons of the beam-bulk model move ahead of the streamer where they ionize the neutral gas. The effect is to increase the streamer velocity relative to the pure drift-diffusion model. The electrons that are emitted ahead of the streamer tip mimic the effect of photoionization as both processes create extra ionization in front of the streamer and increase its velocity. However, the (scaled) electron density produced by the highenergy electrons does not depend on altitude in contrast to photoionization which is affected by the quenching of excited states of the neutral gas.
We further see that the peak electron density and the peak electric field is reduced relative to the pure drift-diffusion model, as also suggested in [15]. The change in altitude (density) in the beam-bulk model has little influence, as it should, because the same scaling laws apply. The only differences are due to statistical fluctuations ahead of the streamer, where it appears that more energetic electrons are created at higher altitudes. We will return to this question in the following section.
The streamer propagation velocities in the two models are shown as functions of time in figure 3 . Once the streamers have developed, the velocity in the beam-bulk model is almost twice as large as the velocity in the pure drift-diffusion model. The peak electric field in the streamer head, E max , is shown in figure 4. When the streamers are formed, E max reaches ∼ E 7 k in the beam-bulk model and appears almost constant whereas the field in the pure drift-diffusion model continues to grow and reaches values above E 9 k at the end of the simulation. The electric field plays an important role for production of high-energy electrons (e.g. [12,14,16]) and for branching of streamers (e.g. [25,26]). We recall that as a first step, we did not take into account photoionization. It is interesting to note that the results presented in [13,26] show similar relative decrease of electric field when the photoionization level increases between a streamer simulated at ground and at more than 30 km altitude. The results presented in [27] show also a streamer acceleration due to photonization that is significant when the background field is really high (∼3 E k ). We conclude, therefore, that it is important to self-consistently include particle acceleration and streamer dynamics to properly reflect streamer properties and their role in TGF generation.

Electron acceleration and the production of runaway electrons
In this section, we discuss the generation of energetic electrons and consider the three beam-bulk streamers of figure 2. In figure 5 we show the unscaled number of electrons that exceed 1 keV as a function of time (scaled). Figure 6 shows the maximum energy of electrons reached in the three    simulations of figure 2. It is important to note that during the discharge ignition for < t 0.2 ns, we observe on figures 5 and 6 a first burst of runaway electrons due to the high electric field close to the electrode. We have checked that this first burst has no influence on the results obtained during the propagation of the discharge studied in this work. After the streamers have formed, figure 5 shows that the production rate (slope) appears independent of altitude, however we see that the total number of electrons increases by a factor ∼4 for every 10 km altitude. This can be understood when we consider the total number of electrons in a streamer. As noted earlier, the density of electrons in a streamer scales as n gas 2 and a volume element scales as n 1 gas 3 . The number of electrons in a streamer therefore scales as n 1 gas . Since the neutral density  There are therefore two aspects that make it comparatively easier with increasing altitude to generate TGFs; one is that the number of source electrons in a streamer increases with altitude, and the other is that absorption of bremsstrahlung photons decreases. The maximum energy grows at the rate of ∼ − 10 keV ns 1 , independently of altitude almost reaching an energy that corresponds to the potential drop of −25 kV available.
We now look further into the acceleration process by analyzing streamer simulations performed at ground altitude for different voltages applied to the point electrode. = qE F d , where F d is the average frictional force experienced by an electron moving through the gas (e.g. [12]), is shown as the red line. It depends on the electric field magnitude, allowing electrons with lower energies to enter the runaway regime for higher electric fields. We see that the energetic beam ahead of the tip is in the runaway regime and takes its origin where the threshold for runaway is minimum and the electric field is maximum. We note here, that at this location, the total electron density is above − 10 cm 11 3 and − 10 cm 14 3 respectively, which is too high to be simulated by real particles. This justifies the approach of the beam-bulk model over other hybrid approaches, when studying electron acceleration in streamers.
In figure 8 we show the total number of runaway electrons in the simulated streamer, N r , as function of time for 5 applied voltages between −25 kV and −35 kV. Here we define the runaway regime as where the force of the electric field on an electron equals the average frictional force at that energy. The number of runaway electrons increases rapidly with decreasing (negative) applied voltage from about 20 at −25 kV to 10 4 at −35 kV. The initial pulse in N r is caused by the rapid increase and overshoot of the peak electric field as  shown in figure 9. The high electric field generates bursts of energetic electrons and, furthermore, leads to a lower threshold energy for the runaway regime, both effects increasing N r . After the initial overshoot, the field remains almost constant in time which means that the simulation system has reached a steady state. We further note that there appears to be a saturation of E max and N r for applied voltages approaching −35 kV. After the overshoot, E max saturates at ∼E k = 8 and N r at ∼10 4 electrons. We expect, as noted earlier, that the N r scales with n 1 gas and is ∼4 times larger at 10 km altitude.
It is important to note that the simulation time is short compared to the more than 6 ns runaway avalanche time for electric fields below E k . Still, we expected N r to increase with time after the initial overshoot. After all, runaway electrons are supposed to be accelerated and, at the streamer tip, the number of runaway electrons is expected to increase. However, N r slightly decreases with time, suggesting that runaway electrons are lost at a rate that balances the inflow at the streamer tip. The answer to this question can be found when re-examining figure 7. Ahead of the tip, the number of energetic electrons in the runaway regime above the red line decreases with distance ahead of the tip because they diffuse to lower energies and are lost to the runaway population. Figure 10 presents the runaway current on the axis for an applied voltage of −35 kV at 1 ns. It shows clearly that the runaway current decreases as the runaway electrons travel away from the tip: it varies form μ 30 A close to the tip to μ 2 A 5 mm ahead. The conventional runaway condition that we have adopted so far is really for electrons with velocity vectors anti-parallel to E, and can be written as ( ) 0 p p d for a particle of energy ϵ p at the position Z p . It does not take into account the effect of the pitch angle relative to E. If we include now the pitch angle in the runaway condition, as in [10] , we obtain: , where θ p is the angle between the velocity vector and the Z-axis. Figure 11 presents the distribution function of runaway electrons defined by  as a function of the perpendicular ( ⊥ V ) and parallel (V z ) velocity components relative to the electric field for an applied voltage of −35 kV at = t 1 ns. Each of the four sub-panels is plotted for an electron population inside a 5 mm thickness section starting respectively at −6, −7, −8 and −9 mm. On the top of each figure, we have plotted the criterion θ  that now depends also on the angle θ p for the average electric field E m , which equals respectively E 2.02 k , E 1.06 k , E 0.8 k and E 0.7 k in the section. On this figure it is clearly seen that the threshold energy defined by θ  increases rapidly with pitch angle and that an important amount of runaway electrons as defined by  do not satisfy θ  . This problem is further aggravated by the fact that the electric field is decreasing away from the tip, leading to an increase in the runaway threshold. These considerations make therefore difficult the evaluation of the flow of cold runaway electrons needed to justify TGF observations. It is still to be noted that the current density from  [12,14,16] to justify TGFs from low energy electron accelerated by streamer tips using extra RREA amplification outside the streamer tip not present in our model. The lower rate we found can be explained by several factors: the electric field is of lower magnitude inside the streamer tip due to the production of high energy electrons itself; the electric field outside the streamer tip is below E k in our work; more runaway electrons are expected to be lost by diffusion and our streamer radius is 1 order of magnitude smaller.

Conclusions
In this paper we have applied for the first time a 'beam-bulk' model to simulate the runaway production from cold electrons accelerated in streamer tips present in lightning leader region. The method is really well adapted to the physics of the process since it models the bulk of low energy electrons present in the body of the streamer using fluid drift-diffusion equations and the high energy electrons that contain the beam of runaway electrons using particles. It allows to observe the production of runaway electrons at the core of their source production in the electron velocity-space. The results suggest that the effect of high energy electrons on streamer propagation is important on the velocity, peak elect ric field, eventual branching and on the production rate of runaway electrons. The effect is similar to that of photoionization but is independent of altitude. It is shown that the model can handle a runaway rate that is few orders of magnitude below what is used by other authors to justify TGF observations. However in our self-consistent simulation, we have observed that a significant amount of runaway electrons are lost by diffusion along their way outside of the streamer tip. This finding makes difficult the estimation of the runaway number needed to justify TGF observations from the present study. Although the model does not yet contain a realistic description of the lightning field, nor extra RREA amplification outside the streamer tip, confrontation of results given by such a model  with TGF observations expected from future space missions like TARANIS and ASIM in 2015/16 can give insight in the source mechanism and in term lead to a better understanding of atmospheric discharges and their impact on our environment [39,40].