Self Consistent Modeling of Relativistic Runaway Electron Beams Giving Rise to Terrestrial Gamma‐Rays Flashes

Terrestrial Gamma Ray Flashes (TGFs) are short bursts of gamma rays occurring during thunderstorms. They are believed to be produced by Relativistic Runaway Electron Avalanches (RREAs). In this paper, we present a new numerical model based on the Particle‐In‐Cell (PIC) method to simulate the interactions between the electromagnetic fields and the electron avalanche self‐consistently. The code uses a cylindrical Yee lattice to numerically solve the electromagnetic fields, a Monte Carlo approach to simulate collisions with air molecules, and a plasma fluid model to calculate the effects of low‐energy electrons and ions. The model is first tested through the reproduction of dispersion relations in a hot plasma. RREAs propagating under a homogeneous background electric field are then simulated. Owing to the self‐consistent nature of this description, we report here new physical properties such as saturation processes in the electron density and in the number of high‐energy electrons, detailed dynamical screening of the electric field in the ion trail of the avalanche, and the associated electric currents. We find that the saturation of RREAs is obtained when the numbers of high‐energy electrons and photons is consistent with those believed to be representative of TGF sources.


Introduction
Terrestrial gamma-ray flashes (TGFs) are short bursts of gamma-rays produced in thunderstorms by avalanches of very energetic electrons called Relativistic Runaway Electrons Avalanches (RREAs).The initiation and sustainability of these energetic electron avalanches is still debated, as they require strong electric fields to accelerate electrons to the runaway regime (∼2.8 • 10 5 V/m, Dwyer et al., 2012), and a substantial number of these energetic electrons to eventually lead to an emission with the intensity of a TGF (∼10 18 photons with high energy ≥1 MeV at the source).Two mechanisms have been considered to produce avalanches conducive to TGFs.One is the Relativistic Feedback Mechanism (Dwyer et al., 2012), that relies on the emission of X-rays and energetic positrons triggering successive avalanches, to sustain the process until the required number of particles is reached.The other is the thermal runaway electron mechanism taking place in streamer discharges and forming RREAs ahead of a leader (Celestin & Pasko, 2011), also called leader-based mechanism, which depends on the propagation of leaders during lightning: the strong electric field and great number of electrons at the tip of a leader could be sufficient to initiate an avalanche large enough to produce a TGF (Celestin et al., 2015).
While many models have simulated accurately the RREA (e.g., see Berge et al., 2022;Dwyer, 2021), and many key features are now known, none have modeled the full interactions between electromagnetic fields and particles.All these radio observations have yet to be reproduced numerically, and full self-consistency was never realized in particle description coupled with electromagnetic fields.
In this article, we describe a new approach of the problem using for the first time in this context a fully hybrid, relativistic, self-consistent Finite Difference Time Domain (FDTD) model, coupling the electromagnetic Particlein-Cell (PIC) method with a Monte Carlo code.In the next part, we detail the structure of the models and methods used, and provide test results to check the ability to model a plasma.In the third part, we describe the simulation results for the case of RREAs obtained using different assumptions on seeding, and report interesting unforeseen features.In the last part, we discuss and interpret the results.

Methods
The main originality of the model is in the coupling between a Monte Carlo method to keep track of collisions and their effects, with a electromagnetic particle-in-cell (PIC) model to obtain a fully relativistic self-consistent description.This section provides a description of the main components of the model.

Monte Carlo Model
In order to keep track of the dynamics of electrons in air in a causal fashion, we use a null-collision-technique Monte Carlo model based on that developed by Celestin and Pasko (2011).The model is 3-dimensional in configuration space, 3-dimensional in velocity space, fully relativistic, and is able to simulate the dynamics of electrons in air over an energy range going from sub-eV to hundreds of MeVs.The air is asummed to be made of 80% nitrogen and 20% oxygen.
Following Moss et al. (2006), 23 excitation collisional processes are taken into account for N 2 (rotational, vibrational, and electronic) and 11 for O 2 including also dissociative processes.The only difference with Moss et al. (2006) is that the dissociation of the oxygen molecule into O( 3 P) and O( 3 S) with a threshold of 14.7 eV is not taken into account in the model.It is expected to play only a minor role.At energies greater than ∼1 keV, which is the energy range of interest in the present study concerning the Monte Carlo model, excitation cross sections are extrapolated logarithmically.Moreover, as angular scattering differential cross sections are not readily available for all excitation processes, the direction of electrons' momenta after excitation collisions is assumed identical to that in the case of elastic scattering (see below).It is worth noting that excitation processes only have an overall minor role in the dynamics of high-energy electrons and are of rather low-probability compared to ionization and elastic collisions.They are described here for the sake of completeness as the model can be used over an energy domain extended to lower energies.
The dynamics of >1 keV electrons in air is dominated by ionization and elastic collisions.Ionization is modeled through the use of the relativistic binary-encounter-Bethe (RBEB) model to obtain an orbital-based analytical singly differential cross sections (DCS) (Kim et al., 2000).It is conveniently integrable analytically and does not need adjustable parameters as it only requires orbital kinetic energies and binding energies of target electrons (Hwang et al., 1996;Santos et al., 2003).The RBEB model provides a seamless coverage over the whole energy domain for primary and secondary electrons well-suited for the description of RREAs (e.g., see Celestin & Pasko, 2010;Celestin et al., 2015;Xu et al., 2015).The direction of the secondary electron's momentum is obtained from the conservation of energy and momentum (e.g., see Celestin & Pasko, 2010).
We use the inverse transform sampling of the DCS to tabulate the energy of secondary electrons as a function of a uniformly distributed random number for all primary electron energies (Moss et al., 2006, Equation 22).However, the ionization DCS shows strong dynamics over the broad energy range of interest and a linear sampling of secondary energies would require an enormous number of points to obtain a sufficient resolution.To avoid it, we empirically found that the following non-linear sampling of the DCS integral was well-suited: where dσ(ε p ,ε s ) dε s dε s is the ionization DCS, σ i (ɛ p ) is the total ionization cross section, and ϵ is set to a value of 10 10 to avoid the argument of the logarithm to reach exactly one.It is now J DCS (ɛ s ) that we sample linearly over N R values (in the present case we use N R = 10, 000, which has proven to be sufficient).Using a random number R ε s uniformly distributed between 0 and 1, the energy index of the secondary electron is then found by taking the integer part of: Although ionization collisions drive the electrons' linear energy loss, elastic collisions significantly impact their dynamics in phase space (e.g., Dwyer, 2010).Below 500 eV, we use differential cross sections obtained  (Kambara & Kuchitsu, 1972;Shyn et al., 1972).The differential cross section for O 2 is assumed identical (for more details on the comparison of angular scattering between N 2 and O 2 at low energy, see Moss et al. (2006), Figure 3, and references therein).
For high-energy electrons (>500 eV) we consider that elastic scattering from molecules is similar to elastic scattering by atoms.Assuming that the electric potential of the nucleus and the atomic electrons result in a screened Coulomb potential, one can show that the angular differential cross section has the following form (Carron, 2006;Dwyer, 2007): where χ is the scattering angle, dΩ is the corresponding differential solid angle, Z is the atomic number (i.e., Z = 7 for nitrogen and Z = 8 for oxygen), r e is the classical electron radius, p is the electron momentum, and a is the atomic screening radius.For the sake of comparison with previous works, we use a ≃ 1.3413Z 1/3 a 0 , where a 0 is the Bohr radius, which is the same value as that used in Dwyer (2007).
In Monte Carlo simulations, the scattering angle χ after one collision is calculated through a uniformly distributed random number R χ between 0 and 1, by finding χ so that: where σ a e (ε) = 2π∫ π 0 dσ e dΩ sin χ dχ is the total elastic cross section for nitrogen or oxygen atoms.Equation 5 may be solved numerically by pre-tabulating R χ .However, high-energy electrons being much more forward-scattered than low-energy electrons, the corresponding sampling would have to be performed cautiously.Instead, it is possible to obtain an accurate solution of Equation 5 through analytical considerations and this is the method we use.
2 and B = ℏ 2 4p 2 a 2 , plugging Equation 4 in 5, one can integrate analytically Equation 5.One finds: where X = 2B cos χ + 1.To solve directly this equation for X, and then for χ, one could use the inverse function of f(W) = W exp(W), where W is named the Lambert W-function.However, the Lambert W-function is uncommon and its use is not necessary in this context as one can make the following approximations.Over the energy range of interest, the term varying as ∼(Bβ 2 + 1)/X in Equation 6 is always dominant over the term in ∼β 2 ln(X ), and Bβ 2 ≪ 1.Hence, one gets: The scattering angle χ can be calculated from R χ through the following formula: Note that this method, although very accurate, is an approximation and the domain of validity for 0 < R χ < 1 is not fully covered.The following condition has to be added: if the argument of the arccosine on the right-hand side of Equation 8 is found to be lower than 1, then χ is set to exactly π.The highest probability of occurrence of the latter condition is reached for electrons with energy ∼3 keV, and corresponds to only 1.33 times over 100,000 elastic collisions on average.For electrons with energies 1 and 10 MeV, the occurrence rates are 1.55 over 1 million and 3.2 times over one hundred million, respectively.
One can calculate the total elastic scattering cross section through the integration of Equation 4: One finds: which is the elastic scattering cross section for nitrogen or oxygen atoms.When calculating the probability of elastic collisions of high-energy electrons with nitrogen or oxygen molecules the molecular cross section should be used: σ m e (ε) ≃ 2σ a e (ε).In addition, for the calculation of the total elastic scattering cross section using Equation 10, from which the collision probability in the Monte Carlo code is found, we use the classical screening radius in the Thomas-Fermi model: a ≃ 0.885Z 1/3 a 0 .Indeed, comparisons of the total elastic scattering cross section given by Equation 10 with the Evaluated Electron Data Library (EEDL) (Cullen et al., 1991) at different energies for nitrogen atoms and with (Itikawa, 2006, Table 3) at 1 keV for nitrogen molecules (using σ m e (ε) ≃ 2σ a e (ε)) have showed a very good agreement for this magnitude of the screening radius in the elastic cross section.

Electromagnetic PIC Model
One of the difficulties of simulating plasmas is to correctly describe the electromagnetic interactions between charged particles.The simplistic solution would be to compute at each timestep all the Coulombian interactions produced by each particle on all the other particles (i.e., solving the N-body problem), but this quickly becomes unbearable when the number of particles becomes relatively high.The PIC method proposes another solution: particles move freely according to the relativistic equations of motions, but the electromagnetic field they create is described and updated over a grid, and is then recomputed onto the particle location via an interpolation technique.The computation gains come from the fact that the number of grid points can be strongly limited and each particle handled in simulations actually represent a great number of real particles (electrons in the present work).This is made possible through a careful discretization of Vlasov equation into characteristic curves resulting in the trajectories of these "macro-particles" or "computer-particles" in phase space.The use of electromagnetic particle-in-cell code to model and simulate plasmas is described in Birdsall and Langdon (1991) and Hockney and Eastwood (1966) (see also the remarkably clear description given by Lehe ( 2014)) The method relies on solving Maxwell's equations, that will be used to update the electric and magnetic fields at each time step.Our model is a cylindrical axisymmetric implementation of the PIC method.We identify four steps in a cycle of the PIC code: • Interpolation of fields onto the particles (Section 2.2.3) • Update of the particle positions using the relativistic equation of motion (Section 2.2.1) • Assignment of the particles charges and currents on the grid (Section 2.2.3) • Update of the electromagnetic fields using Maxwell's equations (Section 2.2.2) We detail these different steps in the following sections.

Equation of Motion of the Electrons
Describing the particle movements is straightforward, as we simply use the equations of motion, in the following form: where E is the electric field, B the magnetic field, q e the particle charge, v i , p i and x i are respectively the particle i velocity, momentum, and position vectors, m e the particle mass, and γ i the Lorentz factor.The use of relativistic equations is required here, as the typical electron speed in a relativistic avalanche is around 0.9c, with c the speed of light in vacuum.To discretize the equations, we use the relativistic Boris algorithm, as the naive approach results in breaking energy conservation and introduces significant errors on the cyclotron frequency.We first define η i = q e |B n+1 |Δt/(2γ i m e ).The Boris algorithm computes the momentum as follows: where xyz is either x, y, or z, i is the index of the electron, n is the time index, E i , B i and F rad,i are respectively the interpolated electric and magnetic field at the location of particle i and the continuous radiative friction applied on this particle.In the present model, the production of individual photons does not affect the dynamics of the particles and this is why a radiative friction is included (especially for electrons with energies >10 MeV).It is taken from Berger et al. (2005).We then update the velocity and position:

Computation of the Electromagnetic Field
The electromagnetic field is computed using a finite-difference time-domain (FDTD) scheme on a grid at each timestep using the following Maxwell-Faraday and Maxwell-Ampere equations: where c is the speed of light in vacuum, μ 0 the permeability of free space, and J is the conduction current density.This current density, as well as the particles densities and charge density are all computed on the same numerical grid as that used for the electromagnetic field.The current density J is calculated through summation of the current assigned onto the grid from computer particles (see Section 2.2.3) and the current densities of ions and low-energy electrons calculated using fluid equations (see Section 2.3).
To initialize the field, we use the Maxwell-Gauss equations: Journal of Geophysical Research: Space Physics with ɛ 0 the permittivity of free space and ρ the charge density.In practice, the electric field is initialized with a Poisson equation, and a homogeneous field is added in the cases presented in this article.The magnetic field is initialized as equal to 0 on the whole domain.Equations 22 and 23 are considered as initial conditions if charge is conserved.If they are true at instant t = 0, they remain true at all instants, provided that the continuity equation for the charge is verified: It is hence critical that the numerical model ensures that the charge of the system is preserved (see Section 2.2.3), otherwise a Poisson equation corrector has to be implemented to establish consistency between the electric field and the charge density (Birdsall & Langdon, 1991, Section 15.6).
The grid used is a Yee lattice, meaning a combination of two grids, one for the cells themselves, or cell centers, and another for the interfaces.The grids are staggered from one another, and thus two sets of coordinates are used.Different properties represented on the grid are placed at different positions relative to the cell centers, which allows for consistency when, for example, updating the magnetic field at a point requires to know the rotational of the electric field at this point.The form of the grid used, as well as the different field components, are shown in Figure 1.
Because of the axisymmetry of the model, we can restrain the computation of the different fields to the derivation of the components represented in Figure 1.If we assume the distance between two cell-centers is equal to Δz (resp.Δr), the distance between a cell center and its interfaces is Δz 2 (resp.Δr 2 ).On this grid, the discretization for the electric field then goes as follows: where E z and E r are the z and r components of the electric field, J z and J r the z and r components of current density, B θ is the θ component of the magnetic field, RG and RIG are arrays with the r coordinates of respectively cellcenters and interfaces.b1 is an intermediate variable, and is not used elsewhere.
The magnetic field is also staggered in time compared to the electric field, which must be taken into account in the Maxwell equations.The general form of the magnetic field equation excluding boundaries is: It should be noted that the staggering in time must be taken into account in the particle movement as the particles move over integer timesteps.To do that, we average B θ n 1/ 2 and B θ n+1/ 2 to have the magnetic field at time nΔt.

Assignment of Charges and Currents and Field Interpolation
The passage from the particles to the cells of the grid is the core of the particle-in-cell method, and must be handled carefully, especially regarding charge conservation.
For the charge assignment, a first order scheme called cloud-in-cell (CIC) is used.The scheme is widespread in PIC-type simulations and described in Birdsall and Langdon (1991) and Hockney and Eastwood (1966).On a 2-D grid, the idea is to distribute the charge among the four nearest cells by multiplying it with what is called a shape factor: each cell has its factor, which is comprised between 0 and 1, and the sum of the four factors is equal to 1. Graphically, the particle can be represented on a 2-D grid as a uniformly distributed "cloud," that is separated into four "sub-charges," as shown in Figure 2.
If we name the cells (i, j), (i + 1, j), (i, j + 1), and (i + 1, j + 1) respectively 1, 2, 3, and 4, we can describe each shape factor by the following formula: where k is the cell index equal to 1, 2, 3 or 4, S k is the shape (or form) factor, Δz and Δr are the interval between cells respectively in the z and r direction, as shown in Figure 2, and z loc k and r loc k are respectively the local coordinates of the particle relative to the corresponding cell.We can write them as such: r loc 1,2 = r j + Δr r p = r j+1 r p (32) where (z p , r p ) is the absolute position of the particle in the simulation domain, and (z i , r j ) are the coordinates of the corresponding cell.
The contribution of particles on the charge density assigned to each cell is the shape factor multiplied by a "local" charge density n k = q p V cellk : where q p is the particle charge and V cellk the volume of a cell.In a 2-D Cartesian grid, V cellk is usually constant.However, in the case of our cylindrical axisymmetrical grid, each cell is a torus of revolution with a square section ΔzΔr and a width depending on the position of the cell on the r-axis.If 0 ≤ r j ≤ Δr 2 , then the shape of the cell is a cylinder.To take into account this geometry when computing the charge density, we use a different volume for each coordinate j on the r-axis: As said previously, the same scheme will be used to interpolate the field on the particle.Note that when interpolating the magnetic field on the particles, we take the average of B n 1/2 and B n+1/2 to obtain the value B n at the current timestep nΔt before making the interpolation.The fields applied to the particle will be the sum of the fields at each of the four nearest cells multiplied by the corresponding shape factors so as to avoid the introduction of spurious self-forces.Indeed, the PIC methods separates the calculation of particle motion and the evolution of fields.The fields are therefore not necessarily physically consistent with the position and motion of particles and spurious behaviors such as a particle exerting a force on itself, the so-called self-force, as well as other inconsistencies can occur.To avoid such problems, the shape factor used for charge and current density assignment to the grid simply needs to be the same as that used for the field interpolation (Hockney & Eastwood, 1966).It can be proven that in this case, self-forces and other inconsistencies cancel out.For this reason, the search for typical inconsistencies is a convenient way to sanity-check PIC simulation codes.
However, we can demonstrate that a current assigned though the CIC scheme does not automatically verify the continuity equation (e.g., see Lehe, 2014).Let us consider a 1D-grid with a cell size of Δz on which a particle of charge q p is moving.The timestep is Δt.We use the cloud-in-cell assignment scheme for both the charge and the current density.The shape factor S(z z p ) is therefore expressed as 1 Δz for |z z p | < Δz and 0 otherwise, with z the position of the point being assigned, and z p the particle position.The situation from t = nΔt to t = (n + 1)Δt is represented in Figure 3.
To validate our scheme, we need to verify that the charge is conserved over all cells.Particularly, we can express charge conservation on the cell-center of position kΔz as follows: We can express ρ and J z : ρ We can write the two terms of the continuity equation: As we can see when we compare Equations 40 and 41, the sum of the two is not equal to 0 in all cases.Therefore, using the CIC scheme for the current assignment is not consistent with charge conservation.
A way to correct this discrepancy is to use a Poisson corrector on the electric field, canceling the difference introduced by current assignment (Birdsall & Langdon, 1991).However, the time cost makes it an inconvenient solution, so instead we use a specific current assignment scheme, described in Villasenor and Buneman (1992).Like in the CIC scheme, the particle is described as a cloud with the shape of a cell and an uniform charge, however the shape factor is not computed using the distance of the particle to the computing points, but rather the amount of charge going through each interface, determined by the volume of the cloud going through it.
Provided that the CFL condition is verified, the movement over one timestep can only affect four interfaces at once (movement relative to the nearest interface point), seven interfaces at once (movement from one interface point to another in either the axial or radial direction), or in the rarest cases ten interfaces at once (movement from one interface point to another in a diagonal direction).
We can quickly check that the CFL condition is verified: in most cases presented in this article, the cell dimension is Δz = Δr = 8 m, and the timestep used is Δt = 3.4 × 10 13 s.The Courant number is computed as such on our grid: with U z and U r the z and r components of the velocity.We can maximize the Courant number by taking c as the value of each component of the speed, thus giving a value of C = 2.5493 × 10 5 , which is way below 1, hence verifying the CFL condition.
The different types of movements are shown in Figure 4. We have verified that the last kind of movement (Figure 4c) is rare enough so that we can neglect them without a significant loss of charge conservation.The various runs performed have shown close to no such movement, with sometimes only one at the very beginning, and none afterward, validating this assumption.

Remapping of Particles
Keeping track of the exponential multiplication of runaway electrons requires to limit the number of superparticles in the simulation.Various techniques have been described in the literature to achieve this goal (Schmalzried et al., 2022).In the present work, we use a simple resampling technique that we refer to as remapping.
Each super-particle amounts for a certain number of electrons, a property called weight and noted W. When a super-particle of weight W ionizes a molecule, the knocked-off electron is set with the same weight as the incident particle, and the event can be viewed as "W electrons have ionized a molecule." When the number of super-particles reaches a threshold number (equal to 10,000 in our simulations), a remapping is performed.It consists in the merging of pairs of particles with one another, resulting in a particle with a weight equal to the sum of the two merged particles.Of course, such manipulation entails significant loss of information, so it cannot be applied simply by randomly merging any two particles.A sorting of particles must be made before the merging, to minimize the loss of information.Depending on the case, different sorting patterns can be used.
In order to save as much information as possible in the energy space, we use a remapping based on neighbors sorted in energy: we remap groups of particles present in the same grid cell, sorting them by energy in this area.It allows us to preserve to a certain extent information on energy and position.

Fluid Model
Keeping track of the dynamics of electrons over the whole energy range of interest would require significant computational resources.In order to shorten simulation runs, all electrons with energies less than 1 keV and all ions are assimilated into a plasma fluid model.The density of low-energy electrons and ions are calculated from particles using form factors discussed previously.The fluid approach for low energy electrons over the space scales considered in this work is justified by the fact that 1 keV electrons slow down and thermalize over length scales much shorter than the grid resolution.For example, the continuous slowing down approximation (CSDA) range of a 10 keV electron in ground-level air is only a few millimeters (calculated from the NIST Standard Reference Database 124, https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html), while the numerical grid used for the fluid model is the same as that used for the PIC code (see previously) that has a resolution on the order of 1 m.
Secondary electrons are explicitly taken into account in the Monte Carlo code, and their contribution is therefore directly accounted for in the fluid electron density.If produced with energies above the ionization potential of molecules, these secondary electrons could results in more ionization.Thus, they are also considered to contribute to the fluid electron density as a linear function of their energy.Indeed, it is well-known that ionizing electrons loose an average energy of ∼34 eV per electron-ion pairs produced in air (e.g., Knoll, 2000, Table 5.1, p. 132), which is approximately constant as a function of energy.Simulations performed with the Monte Carlo code presented here as part of a preliminary work confirm this point.Each electron produced in the Monte Carlo code with energy E lower than 1 keV is thus considered to contribute to an amount of E/(34 eV) electrons and positive ions to the plasma fluid densities.
As the medium is highly collisional, the dynamics of electron and ion densities is considered to evolve according to drift-diffusion equations, as usual in discharge physics (e.g., Bourdon et al., 2007) and also used in the context of RREAs to simulate the effects of low-energy electrons and ions (e.g., Berge et al., 2022;Liu & Dwyer, 2013): where n e , n n , and n p are the electron, negative ion, and positive ion densities, respectively; ν i is the electron-impact ionization frequency, and ν a is the electron attachment frequency considering dissociative attachment producing a negative oxygen atom O (twobody attachment) and a three-body attachment producing a negative oxygen molecular ion O 2 ; and v e , v n , and v p are the electron, negative ion, and positive ion drift velocities, respectively.Because of the magnitude of the densities and the timescales of interest in the present work, we neglect ion-ion recombination, electron-ion recombination, and the diffusion of ions (see Berge et al., 2022, Section 2.1).Given the electron densities obtained in the present study, diffusive fluxes of electrons are much lower than drift fluxes over the smallest length scale considered here, which is on the order of 1 m.The diffusion term for electrons is therefore also neglected.
As already noted in previous studies (Berge et al., 2022;Dwyer & Cummer, 2013), the dynamics of low energy electrons and ions is critical to the production of electromagnetic radiation in the LF range.The electron source terms and transport parameters are obtained from Morrow and Lowke (1997) and ion mobilities from Dhali and Williams (1987) (see Berge et al., 2022, Section 2.1).
As it removes free electrons over short time scales and produce negative ions that respond to the field over longer time scales than electrons, the electron attachment processes are of prime importance in the system studied.We use the two-body and three-body attachment rates as described by Morrow and Lowke (1997, Appendix).However, for fields lower than 10 4 V/m at ground-level, the three body attachment frequency given by Morrow and Lowke (1997) reaches non-physical values (e.g., in comparison with Kossyi et al. (1992), Section 2.4).To avoid this issue and for the sake of simplicity, we cap the three-body attachment frequency at 10 8 N 2 / N 2 0 s 1 , where N is the local air density and N 0 is the ground-level density as shown in Figure 5.
The electric current density obtained from the fluid equations is added to the current density of computer particles when calculating the evolution of the electromagnetic field (see Section 2.2.2).

Domain Boundaries
The domain boundaries are, depending on the case, either conducting walls or periodic.The periodic boundary conditions are used in the test cases, when we only use the PIC part.The conducting boundary condition is one of the simplest, that is implemented simply by considering the electric and magnetic fields equal to zero at the boundaries.It brings several issues, as the electromagnetic waves bounce on these boundaries, leading to non-physical effects if Comparison between three-body attachment rates at ground-level as functions of energy obtained from Morrow and Lowke (1997, Appendix) (red line) and (Kossyi et al., 1992, Section 2.4, Equations 45 and 46) (blue line).For the latter, the air temperature is chosen as T = 273 K and the electron temperature T e is calculated through the Einstein relation using Morrow and Lowke (1997) diffusion and mobility coefficients at any given electric field.The black dashed line is the value chosen in the model presented here for electric fields below which the value given by Morrow and Lowke (1997) exceeds it.
we let the simulation run for too long.However it remains negligible as long as the bouncing EM wave has not reached the RREA, which allows us to observe the electron avalanche for relatively long durations.

Validation of the PIC Code Using the Simulation of an Homogeneous Relativistic Plasma
Additionally to validation tests of the FDTD, Monte Carlo, and fluid parts alone, we have tested the coupling of fields and particles performed by the PIC method.For single particles, no self force was observed (Hockney & Eastwood, 1966).For multiple particles, we have defined a homogeneous plasma of relativistic electrons and static ions.
In this simulation, we fill the domain (z = 500 m, r = 150 m here) with electrons, the velocities of which are drawn from a Maxwellian distribution of temperature T = 2.57 • 10 7 K.The domain boundaries are defined as periodic.
In the first run of this case, we put 40,000 particles with a weight of 2 • 10 14 , giving us 8 • 10 18 electrons.We put as many ions to neutralize the plasma.The initial electromagnetic field is considered equal to zero.We then move some of the electrons to form a hole in the electron population in the (z, r) domain, hence creating a charged location that will serve as an initial perturbation.We let the simulation run for about 10 μs, and observed that the velocity distribution remains Maxwellian.
As noted by many authors, waves present a valuable tool to perform verification tests of plasma models (e.g., Palmroth et al., 2018).Therefore, the main purpose of this test is to verify that we obtain the correct dispersion relation.The dispersion relation of electromagnetic waves in a plasma is described as follows: where ω EM (k) is the angular frequency in rad/s for the electromagnetic waves, ω p the plasma frequency, c the speed of light, and k the wave number.However, the use of the Yee scheme with a staggered grid to compute the various components of the electromagnetic field induces changes in the behavior at higher frequencies and wave numbers approaching the Nyquist limit.It is in fact dictated by the relation dispersion of the Yee scheme (Kilian et al., 2017, p. 4.1): To observe it, we apply a 2-D Fourier transform on time and one dimension of space.We show the result of the Fourier transform of the magnetic field in Figure 6.
On this figure, we also show the two analytical dispersion relations described above, superposed on the Fourier transform.We can clearly see the effect of the dispersion relation for electromagnetic waves in Figure 6 with the plasma frequency as a minimum.We notice a loss of consistency with the dispersion relation for electromagnetic waves as k increases, as expected owing to the numerical dispersion induced by the FDTD computation.The fact that the observed dispersion relation exactly fits the expected numerical dispersion of the Yee scheme demonstrates the correct implementation of the code (Kilian et al., 2017;Palmroth et al., 2018).However, care should be taken when interpreting numerical results for higher frequencies (>50 MHz).

Results
In this section, we will present the results of two simulations performed with the complete hybrid code (Monte Carlo and Particle-In-Cell coupled with the fluid code).Those two simulations are designed to be remotely consistent with the production of TGFs in a thundercloud.Although the acceleration distance and initial number of electrons (or injection rate in the second simulation case, see below) are not fully consistent with the conditions in which TGFs occur, the configuration is simplified to capture interesting basic phenomena.The altitude is 12 km and we assume the air density constant over the whole domain.In both simulations, the electric field will initially be homogeneous, oriented along the z-axis, and strong enough to accelerate the electrons and sustain the avalanche.The particles are super-particles, and remapping will occur when their number reaches 10,000.The simulation domain has (z, r) dimensions of 4,200 m × 1,200 m, on a grid of 500 × 150 cells.The first simulation is a case where all particles are injected at once at the initial time.The second simulation instead uses a continuous injection of particles, in which electrons are injected at a constant rate into the domain.

Instantaneous Injection of Runaway Electrons
Here, we describe the simulation of an electron avalanche made with our complete code.
We inject initially 1,000 particles with a weight of 10 8 each.Their initial position is randomly assigned onto a hemisphere of radius of 42 m at z = 0 m and r = 0 m.We set their initial energy at 1 MeV, their speed is computed from the energy and they are assumed to all go initially in the z-direction.As the avalanche propagates, it will be centered on the z-axis.In the present work, the time step is constrained by the highest probability of collision in the model through the null-collision technique, which yields Δt ≃ 3.40 × 10 13 s.
In this simulation, the background electric field is homogeneous.Its value is 1.1828 × 10 5 V/m (which corresponds to ∼5 kV/cm at ground level).
Figure 7 shows the energy distribution after a steady state is reached.The energy spectrum is consistent with the expected energy distribution in a RREA (e.g., Dwyer et al., 2012).Using a fitting procedure we find that the exponential cut in this case is at ∼8 MeV.

Low-Energy Electron Density
The evolution of low-energy electron density (i.e., the electron density handled in the fluid approach) is shown in Figures 8 and 9.
During the first part of the propagation, the avalanche increases in density rapidly, while expanding slowly, following a typical diffusion trend (proportionally to ̅̅ z √ ) (panels a and b of Figures 8 and 9).At a certain point in time, the electron avalanche reaches a critical density of approximately 10 15 m 3 (panel c), from which the electron density stops increasing, and the avalanche front starts expanding rapidly (panel d).We notice that, while the electron density is at its maximum at the front of the avalanche, it quickly starts to decrease in its tail, while the positive and negative ion densities remain stable (in log scale).We have verified it is mainly due to the 3-body attachment process.

Number of High-Energy Electrons
The number of high-energy electrons can be seen in Figure 10.The number of electrons seems to slow its increase after a time.Comparisons with the profile of electron density (Figure 9) show that the change of behavior happens around the same time the avalanche saturates.
And so, we can determine the number of bremsstrahlung photons emitted by these electrons using Equation 4from Celestin et al. (2015): where N γ is the number of photons produced, N e (t) the number of electrons as a function of time, and 〈ν γ 〉(t) the photon production frequency per electron.At the end of the simulation, we obtain for this simulation 1.5877 × 10 19 photons above 1 MeV produced, which is on the order of the expected number of photons in a TGF (e.g., Celestin et al., 2015, Table 1).

Electric Field
The evolution of the electric field on the z-axis (in the center of the avalanche) is shown in Figure 11.This figure shows the values of the different components of the electric field, as well as the position of the energetic electrons (>1 keV), which gives us an approximation of the position of the avalanche front.Once the avalanche has reached saturation, we observe that the electron density is sufficient to quickly screen the electric field out at the front (panels c and d).However, the depletion of the electron population due to the attachment cancels partially this effect, allowing the remaining electric field to still accelerate the electrons behind the front, which explains the groups of electrons forming there.We also observe the formation of a peak in the electric field at the position where the avalanche saturation was reached (around 1,800 m), negative for E z .In the case of a continuous injection, it would result on an acceleration in the z-direction (see below).

Electric Current
The current produced by the avalanche can be separated into two components: • The conduction current, which corresponds to the current produced by the movement of electrons and other charged particles in time.It is computed as part of the current assignment procedure using the Villasenor scheme (see Section 2.2.3).• The displacement current is produced by the variation of the electric field in time, and computed using the following formula: Journal of Geophysical Research: Space Physics 10.1029/2023JA032278 where j D is the displacement current density, E n and E n 1 are the electric field at time t n and t n 1 , Δt the timestep and ɛ 0 the vacuum permittivity.
The current is plotted in Figure 12.We can see that at the front, both current components seems to be roughly equal to one another (one being positive and the other negative, which explains why the total current is of the order of 10 5 A).We notice that the displacement current is modified ahead of the avalanche.It is to be expected, as the  electrons travel at slower speeds than the speed of light, while electromagnetic waves travel at c.The fluctuation induced here is however too weak to have a visible effect on the electric field plot.
The rise in displacement current appearing from the left of the domain occurs after the rebound of the initial electromagnetic wave on the radial boundary and is believed to be a numerical effect.
Behind the avalanche front, attachment reduces the electron population, so the conduction current starts decreasing (panel c and d).In this simplified simulation, the total current reaches values above 10 5 and even up to 10 6 A. Using a different model, Berge et al. (2022) showed that slow LF pulses reported by Pu et al. (2019) would correspond to peak current of ∼100 kA.

Radiated Field Observed at a Distance
We can compute the magnetic field at a horizontal distance R = 150 km from the point of emission, considering an antenna of height H, using Equation 4 from Uman et al. (1975): where i(z, t R/c) is the conduction current.Here, θ is considered to be equal to π/2 (hence sin θ = 1).
The two integrals of the equation are plotted in Figure 13, and labeled respectively the magnetostatic field and radiation field.When calculating the magnetic field at a distance from the electric current using Uman et al.'s equation, Berge et al. (2022) used the total current.Upon further inspection, only the conduction current should be included in Uman et al.'s integrals as per use of Lorenz retarded potentials (see also Shao (2016)).We can estimate that the error caused by this oversight on the magnitude of the magnetic field reported by Berge et al. (2022) is ∼20%.
Since the ambient electric field is strong, homogeneous and infinite in the simulation, and that initial particles are injected all at once, the form of the magnetic field is not expected to resemble that associated with a TGF.We also note that the value reached by the magnetic field is far too high compared to measurements of slow LF pulses (e.g., Pu et al., 2019), as consistent with the current values presented in the previous section.

Concluding Remarks
This first simplified simulation allows to isolate interesting effects, like the saturation of the densities and the screening of the electric field, as well as a saturation in the increase of the number of runaway electrons.However the configuration remains only remotely reminiscent of the context giving rise to a TGF, because of the homogeneous intense electric field and the instantaneous injection of all the particles at the beginning of the simulation.However, the effects observed will keep their relevance when particles are continuously injected into the domain (next section), hence being affected by the changes induced in the wake of the avalanche.

Continuous Injection of Runaway Electrons
In this section, we present the results of a simulation working with a continuous injection: one particle of energy 1 MeV will be injected in the starting area every 1,000 timesteps.Each particle has a weight of 10 5 electrons.The injection rate is therefore about ∼3 × 10 5 electrons per nanosecond.As in Section 3.1, the background electric field is 1.1828 × 10 5 V/m (∼5 kV/cm at ground level), and homogeneous over the whole the domain.The simulation initially begins with 10 particles, representing 10 6 electrons.The others parameters of the simulations are the same as in the previous part.

Density
The electron density field is shown in Figure 14 (see also the video provided as Movies S1 and S2).We observe here the same behavior as in the simulation with instantaneous injection after the start for the avalanche front.The behavior of the electrons behind the front can now be observed.It seems that, due to the screening of the electric field in the avalanche trail, the electrons cannot move forward anymore.Instead, they are agglomerating near the axis.
The depletion in electron density that was forming behind the avalanche front at the location where saturation is reached (Section 3.1) does not appear in this case (Figure 15) (here around 3,000 m); instead the density starts increasing beyond the previously observed limit, and starts "propagating" backward on the trail of the avalanche.

Electric Field
The electric field is shown in Figure 16.The screening occurs as we observed in the previous simulation: the separation between the avalanche front and its tail becomes apparent in panel d, once the saturation is reached.
Figure 13.Magnetic field computed from the current, using Equation 7from Uman et al. (1975).The radiation and magnetostatic parts of the magnetic field shown here correspond to the two integrals in that equation.The computation is made for a distance from the source of 150 km.The presence of a peak in the negative of the z-component of the field means that the electrons are accelerated and gain energy in this moving structure.This peak increases up to a maximum amplitude as the structure propagates backwards and screens the electric field in its trail.

Number of High-Energy Electrons
We show the number of high-energy (>1 MeV) electrons in Figure 17.It is interesting to see that, even with a continuous injection, this number also reaches a saturation.Calculating the number of bremsstrahlung photons  Journal of Geophysical Research: Space Physics 10.1029/2023JA032278 created in this simulation (following the method in Section 3.1.2),we obtain 1.2940 × 10 19 photons produced after approximately 34 μs, which is again in the range of what we could expect for a TGF.

Concluding Remarks
Self-consistent features such as saturation of the avalanche associated with a collapse of the electric field observed in the case of an instantaneous injection of the primary electrons (Section 3.1) are also observed with a continuous injection.Another feature is observed in the electric field and in the ion and electron densities in the form of a structure moving backward (with respect to the electrons), which can accelerate high-energy electrons.

Discussion
We can see that the energy distribution shown in Figure 7 is close to that reported in previous simulations without taking self-consistent effects into account (e.g., Dwyer et al., 2012).A similar distribution is obtained in the continuous injection case.
For the sake of simplicity, we do not consider the production of photons and positrons explicitly in the present work.Phenomena associated with relativistic feedback (Dwyer, 2003) are therefore not included.Under the electric field applied here (5 kV/cm at ground-level), relativistic feedback should have significant effects beyond a distance of ∼1,200 m (at 12 km) (Dwyer, 2003;Skeltved et al., 2014).However, for a closed system like in our configuration, it seems that this electric field magnitude would correspond to a distance of ∼4,000 m (Pasko et al., 2023), which is close to the size of our simulation domain.Additionally, it is likely that the effects of relativistic feedback would not be significant over the relatively short timescales used in this study.It is also worth noting that the effect of a continuous injection of primary electrons are not unlike the effects of the production of new avalanches by photons and positrons.Similarly to the relativistic feedback, the continuous injection introduces another timescale to the system, additionally to that associated with RREA processes.
In both configurations studied, the saturation of the electron density can be explained by the increase of conductivity associated with a RREA.In the beginning, the increase in the electron and ion densities is controlled by the runaway electron avalanche and characterized by the avalanche time, while the relaxation of the electric field is related to the conductivity in the medium through the Maxwell time.The latter also characterizes the time scale of the defocusing of a charged particle beam through Couloumbian repulsive forces.When the conductivity driven by the low-energy electron density reaches a magnitude such that the Maxwell time approaches and exceeds the avalanche time, the field collapses faster than the time it takes for an avalanche to increase the density.
These two times are defined as followed: where λ is the characteristic avalanche e-folding length, V e is the runaway electrons speed, q e is the electron charge, μ e is the electron mobility, and n e the electron density.Equaling those two timescales, we can thus find the density at which the Maxwell time reaches the avalanche time: ⇒ n e = ε 0 q e μ e τ RREA (53) The term τ RREA depends on the characteristic avalanche length λ, which itself depends on the electric field.For example, it can be found by using Equation 2.3 in Dwyer et al. (2012).We can also use our simulation results by plotting the theoretical increase in density for a given λ and compare it to the increase in the ion density in our simulation (before saturation) or to the peak of the electron density at different moments of time, as the electron Journal of Geophysical Research: Space Physics 10.1029/2023JA032278 density has a faster dynamics due to the attachment processes.The electron density and the theoretical density are shown in Figure 18.For this simulation, we obtain a good fit for λ ≃ 165 m.The speed is assumed constant at V e = 0.89c (Coleman & Dwyer, 2006), with c the speed of light, so τ RREA = 3.779 × 10 8 s.We also have μ e = 0.2719 m 2 /(V.s)Morrow and Lowke (1997); q e = 1.602 × 10 19 C, and ɛ 0 = 8.854 × 10 12 F/m.It yields a saturation density nth = 3.2806 × 10 14 m 3 .
As shown in Figures 18 and 19, the observed saturation density is in good agreement with our computed value.
Following Equation 53 the low-energy electron saturation density can be calculated for various electric fields.The result is shown in Figure 20.
In the cases studied in the present paper, the low-energy electron saturation density is reached at ∼3 × 10 14 m 3 while the runaway electron density is about four orders of magnitude lower, reaching n r,s ∼ 10 10 m 3 at 12 km.It is expected that this runaway electron saturation density would follow a similar trend as the low-energy electron saturation density shown in Figure 20, and their order of magnitude (at 12 km) would hence not change significantly in realistic fields present in streamer zones of negative leaders.
In the conditions of the simulations presented here, the saturation of the RREA occurs when the number of runaway electrons is close to that expected in TGFs (∼10 18 ) and the growth of this number strongly slows down after that point (see Figures 10 and 17).This total number of runaway electrons is related to the radius of the avalanche when saturation is reached.Considering that saturation is reached at the source of TGFs, resulting in the collapse of the electric field, and assuming the number of runaway electrons involved in the TGF production as 10 18 , the size of the avalanche should be on the order of R ∼ ( 3 4π 10 18 n r,s ) ∼ 500 m at 12 km.
We note that the balance between the avalanche time and the Maxwell time leading to a saturation of the electron density is reminiscent of streamer discharges in which a stable balance between the Maxwell time and the ionization time is present in the streamer head for similar reasons (e.g., Wang & Kunhardt, 1990).This saturation mechanism is different from that reported by Luque (2014) using a 1-D electrostatic model of a RREA.In particular, here the saturation in produced by the increase of the conductivity caused by the strong production of low-energy electrons by the RREA, and not by the effect of the runaway electrons on the field.In our simulations, the latter indeed carry much less charge density than thermal electrons.
In the case of a continuous injection, the backward-propagating structure after the front explosion also follows an electric pulse, saturates the electron density in its wake, and screens the electric field behind it.The structure resembles that obtained by Liu and Dwyer (2013) and named relativistic feedback streamer.It is unclear the case presented by Liu and Dwyer (2013) also reached saturation.They indeed obtained ion densities on the same order of magnitude as that reported here, but the runaway electron density seems to be lower.
Those results could also be important in the context of lightning propagation and corresponding radio emissions, as RREAs seem to be able to produce significant densities in the vicinity of lightning discharges.For example, Bourdon et al. (2010) have shown that streamer discharge propagation is greatly influenced by the background electron density.

Conclusions
We have developed a RREA model based on the hybridation of an electromagnetic Particle-In-Cell method and a fluid model in order to simulate selfconsistently the dynamics of relativistic runaway electron beams in air in the context of TGF sources.
The first presented simulation of RREAs are obtained by initially injecting all the particles at once.It instantiates the existence of a threshold low-energy electron density at which the avalanche density almost stops growing and the avalanche radius starts increasing rapidly.This saturation density is found to result from the equality between the characteristic avalanche time and the Maxwell relaxation time, similarly to screening processes occurring in streamer discharges over much smaller length scales.In other words, the conductivity becomes too high for the ambient field to be sustained over a sufficient duration for RREAs to develop.
Interestingly, upon saturation of the RREA, the increase in the number of runaway electrons and associated bremsstrahlung photons slows down  significantly.For the typical radii of the avalanches obtained in our simulations, which are presumably on the same order of magnitude as those possibly existing in the context of lightning discharges, the total number of runaway electrons and photons obtained is surprisingly close to that expected in TGFs.
The electric current also saturates, as well as the radiated magnetic field.The estimated remote radiated fields 150 km from the source are more than ten times stronger than in the observations of slow LF pulses.However, it should be noted that the acceleration region geometry used here was chosen for the sake of simplicity.Comparisons with radio observations of TGFs such as slow LF pulses can bring additional constraints on the geometry and the injection rates.
Simulations with continuous injection have shown a similar saturation with screening of the electric field, as well as the formation of a backward-propagating structure from the position where the avalanche front exploded.Given the densities and fields at play, it is conceivable that the observed structure impacts further propagation of the lightning leader.
More simulations with inhomogeneous electric fields and time-variable electron injections (either ad hoc or controlled by relativistic feedback) are to be done to better reproduce the dynamics of self-consistent avalanche mechanism in a TGF context.

Figure 1 .
Figure1.Representation of a cylindrical Yee lattice.The solid black lines represent the cell-centered grid, while the dashed black lines represent the interface grid.The interface grid is indexed using "half-coordinates" (+1/2).The different components of the computed fields are represented in red at their respective position on the lattice.The charge density ρ is scalar, and computed on cell centers.The components of the current density J r and J z are positioned at the interfaces.

Figure 2 .
Figure 2. Deposition of the charge of one particle onto the grid.The solid lines represent the "cell grid."The dashed lines represent the "interface grid."The square around the particle has the dimensions of one cell.The four areas that form the square are used to attribute the values of the four form factors that are to be used to assign the charge to the four marked cells.

Figure 3 .
Figure 3. Representation of the movement of a particle at position z p on a 1D-grid.The green area represents the shape of the cloud of the particle.On the grid are represented the cell-center positions (k′Δz) on which are computed the charge density ρ n′ k′ and the interface positions ((k′ + 1/2)Δz) on which are computed the current density J z n′+1/ 2 k′+1/ 2 .

Figure 4 .
Figure 4. Three of the different types of movements considered in the Villasenor scheme.The dashed lines are the interfaces (forming the interface grid), the red lines and associated arrows show the interfaces through which a current is produced.The two squares represent the old and current position of the particle cloud.The movements here are exaggerated for clarity.(a) four-boundary movement; (b) seven-boundary movement, on the horizontal direction; (c) ten-boundary movement.The two interfaces that are affected without touching the clouds in this movement are so because of the decomposition of the movement.

Figure 5 .
Figure5.Comparison between three-body attachment rates at ground-level as functions of energy obtained fromMorrow and Lowke (1997, Appendix)   (red line) and(Kossyi et al., 1992, Section 2.4, Equations 45 and 46) (blue line).For the latter, the air temperature is chosen as T = 273 K and the electron temperature T e is calculated through the Einstein relation usingMorrow and Lowke (1997) diffusion and mobility coefficients at any given electric field.The black dashed line is the value chosen in the model presented here for electric fields below which the value given byMorrow and Lowke (1997) exceeds it.

Figure 6 .
Figure 6.2-D Fourier transform of the magnetic field.The two dimensions used are time and the z coordinate.To obtain the wave number, we averaged the field for each z coordinate on all r.The Fourier transform has been done with 10,000 files each separated in time by 3.34 × 10 9 s.The domain extends along the z-coordinate over 500 m.The plasma frequency ω p = 2.6840 × 10 7 rad/s is obtained from the electron density used in this case.

Figure 8 .
Figure 8. Cross-sectional view of the low-energy electron density [log 10 (particle/m 3 )] at different times.The high-energy electrons appear as black dots.The time passed since the beginning of the run is indicated at the top of each plot.The r-axis has been mirrored for clarity, as the model is axisymmetric.The high-energy electrons are shown in Cartesian coordinates (x, z), to avoid doubling the number of points.

Figure 7 .
Figure 7. Energy distribution of the electrons above 1 keV, smoothed using a mean over 3.4014 × 10 7 s (i.e., 10 files), for a simulation with initial injection of all the electrons.

Figure 9 .
Figure 9. Electron and ion densities along the z-axis at different times for initial injection (the densities shown in this figure are averaged in the r-dimension between r = 0 m and r = 67.2m).The time passed since the beginning of the run is indicated above each plot.The blue, red, and black lines represent respectively the density of electrons, negative ions, and positive ions.

Figure 11 .
Figure 11.Electric field along the z-axis (center of the avalanche) at different times.The r and z components of the electric field are shown as a function of z, as well as the (x, z) particle coordinates.E z is shown on the axis, E r is shown at r = 8 m.

Figure 10 .
Figure 10.Number of electrons with an energy greater than or equal to 1 MeV present and produced since the beginning of the simulation as a function of time.

Figure 12 .
Figure12.Absolute value of the total electric current and its components (conduction and displacement) produced at different moments of time for the case of initial injection.

Figure 14 .
Figure 14.Cross-sectional view of the low-energy electron density at different moments of time (presented here in decimal logarithm).Simulation with a continuous injection of particles.Decimal logarithm of the density field of electrons and highenergy electron positions (black dots) are shown.The high-energy electrons are shown in Cartesian coordinates (x, z), to avoid doubling the number of points.

Figure 15 .
Figure 15.Electron and ion densities along the z-axis at different times for continuous injection (the densities shown in this figure are averaged in the r-direction between r = 0 m and r = 67.2m).The time passed since the beginning of the run is indicated above each plot.The blue, red, and black line represent respectively the density of electrons, negative ions and positive ions.

Figure 16 .
Figure 16.Electric field near the axis r = 0 m (center of the avalanche) at different times for a simulation with continuous injection.The r and z component of the electric field are shown as a function of z, as well as the (x, z) particle coordinates.E z is shown on the axis, E r is shown at r = 8 m.

Figure 17 .
Figure 17.Number of electron with an energy higher than 1 MeV as a function of time, for the simulation with continuous injection.

Figure 18 .
Figure 18.Electron density relative to position for initial injection (blue), and fitted density profile (red).The dash-dotted line is the saturation density computed with Equation 53.

Figure 20 .
Figure 20.Saturation density as a function of the ambient electric field derived from Equation53.In the axis' legends, N stands for the local air density and N 0 stands for the air density at ground-level.

Figure 19 .
Figure 19.Zoom in density of Figure 18, and time-centered around the moment the saturation is reached.