Nini Berge

Sebastien Celestin
email: sebastien.celestin@cnrs-orleans.fr

Matthieu Garnung

Wei Xu

Robert A Marshall

Steven A Cummer

Matthieu Garnun

Modeling Low-Frequency Radio Emissions From Terrestrial Gamma Ray Flash Sources

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

 [START_REF] Fishman | Discovery of intense gamma-ray flashes of atmospheric origin[END_REF]discovered bright, short bursts of gamma-rays emerging from thunderclouds on Earth in the data from the Burst and Transient Source Experiment (BATSE) instrument aboard the Compton Gamma-Ray Observatory. These sub-millisecond bursts were dubbed Terrestrial Gamma-ray Flashes (TGFs), and have since been observed by several spacecraft, notably RHESSI [START_REF] Smith | Terrestrial gamma-ray flashes observed up to 20 mev[END_REF], Fermi [START_REF] Briggs | First results on terrestrial gamma ray flashes from the Fermi Gamma-ray Burst Monitor[END_REF], AGILE [START_REF] Marisaldi | Detection of terrestrial gamma ray flashes up to 40 MeV by the AGILE satellite[END_REF], and the instrument ASIM mounted on the International Space Station (Østgaard et al., 2019).

TGFs are thought to be bremsstrahlung photons originating from the electrons in relativistic runaway electron avalanches (RREAs), either formed in the large-scale thundercloud electric field (e.g., [START_REF] Dwyer | Source mechanisms of terrestrial gamma-ray flashes[END_REF], or in smaller scale inhomogeneous fields produced by lightning leaders (e.g., [START_REF] Celestin | Variability in fluence and spectrum of high-energy photon bursts produced by lightning leaders[END_REF]. High-energy electrons moving in air can "run away" in high electric fields; that is, above a threshold field and particle energy, the electrons will gain more energy from the field than they lose through interactions with the surrounding air. As these runaway electrons propagate, they will transfer enough energy for surrounding electrons to also run away, and thereby cause an avalanche effect (e.g., [START_REF] Gurevich | Increase of the electric field in head-on collisions between negative and positive streamers[END_REF]. The exact mechanism providing highly energetic electron seeds to initiate TGF-creating RREA processes, as well as exactly how this connects to lightning activity, remains disputed (e.g., [START_REF] Dwyer | High-energy atmospheric physics: Terrestrial gamma-ray flashes and related phenomena[END_REF]. It has however been demonstrated that TGFs are correlated to intra-cloud lightning activity [START_REF] Stanley | A link between terrestrial gamma-ray flashes and intracloud lightning discharges[END_REF], more specifically to the initial stages of upward negative leaders (e.g., [START_REF] Cummer | Lightning leader altitude progression in terrestrial gamma-ray flashes[END_REF][START_REF] Lu | Lightning mapping observation of a terrestrial gamma-ray flash[END_REF][START_REF] Shao | A closer examination of terrestrial gamma-ray flash-related lightning processes[END_REF].

While the process of linking satellite observations of TGFs to their parent lightning relied on radio observations, some specific radio signatures have been linked directly to TGF activity. Most notable are Energetic In-cloud Pulses (EIPs; e.g., [START_REF] Lyu | Ground detection of terrestrial gamma ray flashes from distant radio signals[END_REF] and so-called slow low-frequency (LF) pulses. In this paper, we will focus on the latter. Slow LF pulses were first reported by [START_REF] Cummer | The lightning-tgf relationship on microsecond timescales[END_REF], with evidence suggesting that these are a direct signature of TGF production (see [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF] for a detailed observational study of six cases). This idea was further explored by [START_REF] Dwyer | Radio emissions from terrestrial gamma-ray flashes[END_REF], where they use an analytical approach to calculate the expected radio emissions from TGF sources. A recent work by Tilles et al. (2020) show in one particular observation how emissions presumably associated with TGFs related to one another, in particular EIPs, narrow Abstract Relativistic runaway electron avalanches (RREAs) occur when electrons in electric fields in air reach energies above which they gain more energy from the electric field than they lose to collisions with the surrounding atmosphere. RREAs are known to happen in the electric fields in thunderstorms, and are considered to be the mechanism responsible for producing Terrestrial Gamma-ray Flashes (TGFs). As RREAs propagate, they leave a trail of low-energy electrons and positive and negative ions behind. These populations of charged particles will carry currents as they move in the thunderstorm electric field. In the present work, we model the charged species left behind by the propagating RREA, and the resulting radio emissions in the context of injection of thermal runaway seed electrons by a leader. We find that for certain initial conditions, these radio emissions match the slow low-frequency (LF) pulses that have previously been observed concurrently with TGFs. This confirms that the slow LF pulses are likely generated directly by the TGF source itself, as has been previously suggested using a different TGF production model. Slow LF pulses may therefore potentially be used to infer characteristic properties of TGF sources.

BERGE ET AL.

In the present study, we will use a multifluid plasma model to find the current generated by ions and low energy electrons left behind a propagating RREA from first principles. This model handles the movement of charged species and the basic chemical processes they experience. It is non-relativistic, but allows self-consistent calculation of the electrostatic field, and thereby takes the non-linearity of the system into account. We then compare the results with measurements of slow LF pulses reported by [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF].

Methods

Fluid Model

As a RREA progresses through air, it leaves behind a trail of positive ions and low energy electrons. We model their transport using a modified multifluid plasma model similar to those described in for example, [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF] and Ihaddadene and Celestin (2015), which were originally designed for modeling streamer discharges. The model is based on solving the drift-diffusion equations for the charged species, coupled with Poisson's equation:

𝜕𝜕𝜕𝜕𝑒𝑒 𝜕𝜕𝜕𝜕 + ∇ ⋅ 𝜕𝜕𝑒𝑒 ⃗ 𝑣𝑣𝑒𝑒 -𝐷𝐷𝑒𝑒∇ 2 𝜕𝜕𝑒𝑒 = -𝑆𝑆 - 𝑒𝑒 + 𝑆𝑆𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 (1) 𝜕𝜕𝜕𝜕𝑝𝑝 𝜕𝜕𝜕𝜕 + ∇ ⋅ 𝜕𝜕𝑝𝑝 ⃗ 𝑣𝑣𝑝𝑝 = 𝑆𝑆 + 𝑝𝑝 + 𝑆𝑆𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 (2) 𝜕𝜕𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 + ∇ ⋅ 𝜕𝜕𝜕𝜕 ⃗ 𝑣𝑣𝜕𝜕 = 𝑆𝑆 - 𝑒𝑒 (3) ∇ 2 Φ = - 𝑞𝑞 𝜀𝜀0 (𝑛𝑛𝑝𝑝 -𝑛𝑛𝑛𝑛 -𝑛𝑛𝑒𝑒) (4)
where n, 𝐴𝐴 𝐴 𝐴𝐴 , D e , and Φ are the number density, drift velocity, diffusion coefficient, and electric potential, respectively, while q and ɛ 0 are the elementary charge and the vacuum permittivity. The subscripts e, n, and p represent electrons, negative ions, and positive ions, respectively. S + and S -represent the rates of production by electron impact and electron loss by attachment. The latter consists in two-body attachment (or dissociative attachment: one of the oxygen atoms taking on a negative charge) and three-body attachment processes, which involves a third body and leads to a negatively charged oxygen molecular ion. Under the relatively low field magnitudes considered in this study only three-body attachment is significant and the only negative ions are 𝐴𝐴 O - 2 . Except for S RREA , source terms and transport parameters (mobility and diffusion coefficients) for electrons are calculated based on analytical formula provided by [START_REF] Morrow | Streamer propagation in air[END_REF]. Ion-ion recombination, as well as electron-ion recombination are neglected here due to the relatively low density of ions and the short duration of the event overall. Indeed, the reaction rates for ion-ion and electron-ion recombinations would be on the order of 2 × 10 -7 cm 3 /s (e.g., see [START_REF] Kossyi | Kinetic scheme of the non-equilibrium discharge in nitrogen-oxygen mixtures[END_REF]. In the present work, the positive ion density can reach ∼10 14 m -3 (see Figure 6), leading to an electron-ion recombination characteristic time scale of ∼50 ms. The negative ion density can reach ∼10 15 m -3 , leading to an ion-ion recombination characteristic time scale of ∼5 ms for the evolution of the positive ion density. Both these timescales are much longer than the duration of processes considered in this paper. The mobility of negative and positive ions is taken as μ p = μ n = 0.26/P, where P is the local air pressure expressed in torrs [START_REF] Dhali | Two-dimensional studies of streamers in gases[END_REF]. The term S RREA represents the production rate of electron-ion pairs by runaway electrons (see below). The model is axisymmetric and space variables are reduced to r as the orthogonal distance from the axis and z as the distance along the axis. Additionally, the model is electrostatic, and as such does not include relativistic effects.

Using this model, the densities control the electric field, which impacts the source terms, S, and the further transport of particles. The modeled scenario begins with the injection of thermal seed runaway electrons by a leader into an ambient electric field with a time distribution set as an initial condition. As such, a large number of electrons (around 10 15 , but this parameter is varied in different simulations-note that this number of runaway electrons is consistent with observed TGF fluences when considering leader potential drops of ∼200 MV BERGE ET AL.

10.1029/2021JD036040 3 of 14 [START_REF] Celestin | Variability in fluence and spectrum of high-energy photon bursts produced by lightning leaders[END_REF]) are injected into an electric field of about 1.3 times the RREA threshold field, extending over the first 4 km of a 6 km long simulation domain.

We have added a term for the production rate of electrons and positive ions due to the RREA, S RREA , to the original model. The densities of positive ions, n p , and electrons, n e , created by a single RREA at location z(t) can be expressed as:

𝑛𝑛𝑝𝑝(𝑧𝑧(𝑡𝑡)) ≃ 𝑛𝑛𝑒𝑒(𝑧𝑧(𝑡𝑡)) = 𝛼𝛼𝛼𝛼𝑟𝑟𝑟𝑟𝑛𝑛(𝑡𝑡) 𝜋𝜋 ⋅ 𝑅𝑅 2 𝐷𝐷 (5)
where R D is the avalanche radius, α is the total number of electrons produced per unit length (over all ionization generations) (α ≃ 1,500 m -1 at 12 km) (e.g., [START_REF] Dwyer | Low-energy electron production by relativistic runaway electron avalanches in air[END_REF], and t) is the number of runaway electrons, with N 0 as the initial number of runaway electrons, and ξ(t) the number of avalanche lengths depending on the RREA history.

N run (t) = N 0 ⋅ e ξ(
𝜉𝜉(𝑡𝑡) = ∫ 𝑡𝑡 0 𝜈𝜈𝑟𝑟𝑟𝑟𝑟𝑟 ( 𝑡𝑡 ′ ) 𝑑𝑑𝑡𝑡 ′ = ∫ 𝑡𝑡 0 𝑣𝑣𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 𝜆𝜆 (𝑅𝑅 (𝑧𝑧 (𝑡𝑡 ′ ))) 𝑑𝑑𝑡𝑡 ′ (6)
where ν run is the runaway electron production frequency, v RREA ≃ 0.89 c is the mean propagation speed of the RREA, and λ is the runaway avalanche characteristic length [START_REF] Coleman | Propagation speed of runaway electron avalanches[END_REF].

The radius of the avalanche at a certain time can be estimated by:

𝑅𝑅𝐷𝐷(𝑡𝑡) = √ 𝛽𝛽𝐷𝐷 𝑣𝑣𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 𝑧𝑧(𝑡𝑡) + 𝑅𝑅 2 0 ( 7
)
where D is the perpendicular diffusion coefficient [START_REF] Dwyer | Diffusion of relativistic runaway electrons and implications for lightning initiation[END_REF] and R 0 is the initial avalanche radius. The runaway electrons will be spread radially according to a probability distribution, so we chose one e-folding length from the center to define R D . In doing so, we find the coefficient β to be approximately 1/4 through fitting of Monte Carlo simulation results.

The differential equations governing R D (t) [START_REF] Qin | On the inception of streamers from sprite halo events produced by lightning discharges with positive and negative polarity[END_REF] and ξ(t) can be written as:

𝑑𝑑𝑑𝑑𝐷𝐷(𝑡𝑡) 𝑑𝑑𝑡𝑡 = 𝛽𝛽𝐷𝐷(𝐸𝐸(𝑧𝑧(𝑡𝑡))) 2 ⋅ 1 𝑑𝑑𝐷𝐷(𝑡𝑡) (8) 𝑑𝑑𝑑𝑑(𝑡𝑡) 𝑑𝑑𝑡𝑡 = [ 𝜆𝜆(𝐸𝐸(𝑧𝑧(𝑡𝑡))) 𝑣𝑣𝑅𝑅𝑅𝑅𝐸𝐸𝑅𝑅 ] -1 (9)
with R D (t = 0) = R 0 and ξ(t = 0) = 0. Using the forward Euler method, we can therefore calculate n p (z, t) ≃ n e (z(t)) by advancing

𝐴𝐴 𝐴𝐴 𝑘𝑘+1 𝐷𝐷 = 𝐴𝐴 𝑘𝑘 𝐷𝐷 + Δ𝐴𝐴𝐷𝐷 = 𝐴𝐴 𝑘𝑘 𝐷𝐷 + 𝛽𝛽 𝐷𝐷 2 ⋅ 1 𝐴𝐴 𝑘𝑘 𝐷𝐷 Δ𝑡𝑡 and 𝐴𝐴 𝐴𝐴 𝑘𝑘+1 = 𝐴𝐴 𝑘𝑘 + [ 𝜆𝜆∕𝑣𝑣𝑅𝑅𝑅𝑅𝑅𝑅𝐴𝐴 ] -1 Δ𝑡𝑡 ,
where k is an index representing the time iteration, and Δt is the time step of the simulation.

Rather than looking at a single avalanche as above, we need to take continuous injection of RREA electrons into account. The quantities ξ and R D then become dependent on both z and t, as they both depend on the history of the electric field along z(t). The production rate of electron-ion pairs produced by runaway electrons is modeled by the following source term:

𝑆𝑆𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 = Δ𝑛𝑛𝑒𝑒(𝑧𝑧𝑧 𝑧𝑧) Δ𝑧𝑧 = Δ𝑛𝑛𝑝𝑝(𝑧𝑧𝑧 𝑧𝑧) Δ𝑧𝑧 = 𝐼𝐼 ( 𝑧𝑧 - 𝑧𝑧 𝑣𝑣𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 ) 𝑒𝑒 𝜉𝜉(𝑧𝑧𝑧𝑧𝑧) ⋅ 𝛼𝛼 𝜋𝜋 ⋅ 𝑅𝑅 2 𝐷𝐷 (𝑧𝑧𝑧 𝑧𝑧) ( 10
)
where seed electrons are added over time at the location z = 0, according to the injection function I(t). Since standard RREA parameters are used in this source term, although their effects is only modeled using a fluid approach in the present paper, the energy of injected electrons can be understood as distributed according to a standard RREA spectrum with a ∼7 MeV cutoff.

Throughout this paper, we use a Gaussian injection function, given by

𝐴𝐴 𝐴𝐴(𝑡𝑡) = 𝑎𝑎 ⋅ 𝑒𝑒 -(𝑡𝑡-𝑡𝑡 𝜇𝜇 ) 2 ∕𝑡𝑡 2 𝑖𝑖𝑖𝑖𝑖𝑖 ,
where the duration of the injection, t inj , is a free parameter. Because there is no factor of 1/2 in the exponent, 𝐴𝐴 𝐴𝐴𝑖𝑖𝑖𝑖𝑖𝑖 = √ 2 ⋅ 𝜎𝜎 . An injection time of 15 μs therefore corresponds to a standard deviation σ of 10.6 μs. The mean time, t μ is set as 3 ⋅ t inj .

The normalization factor, a, is set so that the function corresponds to a chosen number of initial electrons. A Gaussian time profile of the injected runaway electron flux is used to have a symmetric bell curve characterized by a small number of parameters. Figure 1 shows how the information is propagated through time, t, and space, z. As the RREA moves at the speed v RREA = 0.89c, the time the RREA uses to cross a simulation grid cell of size Δz is 𝐴𝐴 𝐴𝐴𝐴𝐴 = Δ𝑧𝑧 𝑣𝑣 𝑅𝑅𝑅𝑅𝑅𝑅𝐴𝐴 = Δ𝑧𝑧 0.89𝑐𝑐 . For all grid node locations z j , R D and ξ are transported as:

𝑅𝑅 𝑘𝑘+1 𝐷𝐷 𝑗𝑗 = 𝑅𝑅 𝑘𝑘 𝐷𝐷 𝑗𝑗-1 + 𝛽𝛽𝐷𝐷 𝑘𝑘 𝑗𝑗-1 2 ⋅ 𝑅𝑅 𝑘𝑘 𝐷𝐷 𝑗𝑗-1 ⋅ 𝛿𝛿𝛿𝛿 (11) 𝜉𝜉 𝑘𝑘+1 𝑗𝑗 = 𝜉𝜉 𝑘𝑘 𝑗𝑗-1 + 𝑣𝑣𝑅𝑅𝑅𝑅𝑅𝑅𝑅𝑅 𝜆𝜆 𝑘𝑘 𝑗𝑗-1 ⋅ 𝛿𝛿𝛿𝛿 ( 12
)
where the diffusion radius R D , the diffusion coefficient D, the exponential term ξ, and the runaway avalanche characteristic length λ, each depend on their values at previous time steps, along their path of propagation.

The multifluid plasma modeling based on Equations 1 to 4 makes use of a variable time step Δt. This time step is based on the electrodynamic state of the system, which is governed by low-energy electrons (see e.g., [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF]. We emphasize that in the modified model used in the present work, we also determine Δt on the basis of the electrodynamics of the system, while the RREA parameters R D and ξ (Equations 11 and 12) evolve through the use of a longer time step δt, which is based on the grid size Δz (see above). In practice, Equation 11 sets a lower limit on the initial RREA radius chosen, which must be significantly greater than

𝐴𝐴 √ 𝛽𝛽𝛽𝛽Δ𝑧𝑧 2𝑣𝑣 𝑅𝑅𝑅𝑅𝑅𝑅𝐴𝐴
. For a typical value of Δz ≃ 1 m, this implies R D ≥ 10 m. In the following, we use R 0 = 100 m. As part of a preliminary work, we have verified that the initial value of the radius does not affect the results significantly as long as this value is greater than a few times that lower limit.

The electric current i(z, t) is calculated through the conduction current associated with the motion of charged species modeled by Equations 1-3 and the displacement current obtained through the evolution of the electric field given by Equation 4:

𝑖𝑖(𝑧𝑧𝑧 𝑧𝑧) = ∫ 𝑆𝑆 ( ⃗ 𝑗𝑗 + 𝜀𝜀0 𝜕𝜕 ⃗ 𝐸𝐸 𝜕𝜕𝑧𝑧 ) ⋅ ⃗ 𝑑𝑑𝑆𝑆 ( 13
)
where 𝐴𝐴 ⃗ 𝑗𝑗 = ∑ 𝑠𝑠 𝑞𝑞𝑠𝑠𝑛𝑛𝑠𝑠 ⃗ 𝑣𝑣𝑠𝑠 is the conduction current density (the subscript s labels the charged species taken into account) with 𝐴𝐴 𝐴 𝐴𝐴𝑠𝑠 = 𝜇𝜇𝑠𝑠 𝐴 𝐸𝐸 , 𝐴𝐴 ⃗ 𝐸𝐸 is the electric field, and S the cross-sectional area of the simulation domain.

Radio Emission and Propagation

Analytical Formulation

We propagate the radio signal using the analytical formulation from Uman et al. (1975). This is justified by the mostly axial geometry of the RREA system and the fact the observation distance is much greater than the characteristic size of system. The source of the electric current, i(z, t), is modeled as a vertical antenna placed a certain ground distance D and altitude H away from the receiver (see Figure 2) above a perfectly conducting ground plane. The magnetic field can then be calculated using:

𝐵𝐵𝜙𝜙(𝐷𝐷𝐷 𝐷𝐷) = 𝜇𝜇0 2𝜋𝜋 [ ∫ 𝐻𝐻 0 sin 𝜃𝜃 𝑅𝑅 2 𝑖𝑖(𝑧𝑧𝐷 𝐷𝐷 -𝑅𝑅∕𝑐𝑐)𝑑𝑑𝑧𝑧 + ∫ 𝐻𝐻 0 sin 𝜃𝜃 𝑐𝑐𝑅𝑅 𝜕𝜕𝑖𝑖(𝑧𝑧𝐷 𝐷𝐷 -𝑅𝑅∕𝑐𝑐) 𝜕𝜕𝐷𝐷 𝑑𝑑𝑧𝑧 ] ( 14
)
where R is the distance from the altitude element dz to the receiver (see Figure 2).

The effects of ground attenuation at different distances are estimated by applying first-order Butterworth filters and a scaling factor, depending on the distance. A block diagram of the filters used can be seen in Figure 3. These filters have been designed to approximate empirical values for ground attenuation over land (Radiocommunication Sector of ITU, 2007).

Comparison With an FDTD Model

To validate our approach of using Uman et al. (1975)'s analytical formulation for the propagation of the radio signal from our simulated current (see Section 2.2.1), we have compared our results with a model based on the Finite-Difference Time-Domain (FDTD) method. A detailed description of the FDTD model can be found in [START_REF] Marshall | An improved model of the lightning electromagnetic field interaction with the D-region ionosphere[END_REF] and [START_REF] Marshall | Elve doublets and compact intracloud discharges[END_REF]. An FDTD model is much more comprehensive than the simple analytical solution we use, taking the ground, atmosphere, and ionosphere into account but also the spherical geometry of the Earth. Although more accurate, it is consequently much more computationally demanding. Moreover, the level of accuracy needed is determined by the proof-of-concept focus of the present paper rather than a detailed study from which we could extract the value of physical parameters at the source precisely. For this reason, we only use the FDTD model as a validation for the analytical approach.

In Figure 4, we present the result of the model comparison. The current generated by the fluid model presented in Section 2.1 for the case shown in Figure 7 has been approximated by a Gaussian fit of the current in time and space, given by

𝐴𝐴 𝐴𝐴 = -𝐴𝐴0 ⋅ 𝑒𝑒 -(𝑧𝑧-𝑧𝑧 𝜇𝜇 ) 2 ∕𝑧𝑧 2 𝜎𝜎 ⋅ 𝑒𝑒 -(𝑡𝑡-𝑡𝑡 𝜇𝜇 ) 2 ∕𝑡𝑡 2 𝜎𝜎
, where I 0 = 74.2 kA, z μ = 12 km, z σ = 695 m, t μ = 54 μs, and t σ = 15 μs. Doing this does not change the analytically propagated signal significantly, but eases the input into the FDTD model.

Using the default setup [START_REF] Marshall | An improved model of the lightning electromagnetic field interaction with the D-region ionosphere[END_REF], the FDTD simulations are performed for radio propagation over seawater. The chosen source-to-observation distance is set to 150 km. The analytical result is presented in two versions-one with no ground attenuation (i.e., perfectly conducting ground), and one with filters approximating propagation over land. As expected, the FDTD result falls between the two analytical cases (Figure 4), demonstrating the validity of our approach. A diagram representing the filters applied to our simulation results to approximate ground attenuation (first group of filters, left), and to mimic the detector response of the LF sensor at FT (second group of filters, right). The small gray boxes represent scaling factors applied to the signal. All Butterworth filters are first order, while the Chebyshev filters have the order noted. The single pole 1-MHz Butterworth filter acts as a differentiator, which is part of the detector response. For ground attenuation at distances not mentioned here, the distance is rounded to the closest of these three values for the purpose of ground attenuation application only.

Slow LF Pulse Data and Detector Response

In this study, we compare simulated radio emissions to data of three slow low frequency (LF) pulses occurring simultaneously with TGFs [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF]. The slow LF pulses were collected from an orthogonal pair of LF magnetic field coils placed at the Florida Institute of Technology (FT) in the United States. The sensor has an operating bandwidth of approximately 1-300 kHz [START_REF] Cummer | The lightning-tgf relationship on microsecond timescales[END_REF]. To compare simulation results with events observed by this sensor, we apply a series of Butterworth and Chebyshev filters specifically designed to match the detector response, see Figure 3. A graph of the sensor frequency response can be found in [START_REF] Dwyer | Radio emissions from terrestrial gamma-ray flashes[END_REF].

The three slow LF pulses we use were observed on 3 August 2010, 10 August 2014, and 4 September 2015. They are denoted as 20100803, 20140810, and 20150904, respectively based on their observation dates throughout the paper. For an in-depth discussion of the events, see [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF].

Model Flowchart

Figure 5 gives a schematic overview of the flow of information during the simulations. A source term for electron injection is used as input to the multifluid plasma model (see Section 2.1), which gives the current from the moving charged species as an output. This current serves as the source for the magnetic field propagated using the analytical formulation in Section 2.2.1, giving the radio emission at a chosen distance. Filters that mimic the ground attenuation and the response of the antennas with which we compare results are then applied in order to obtain the final signal.

Results

The following results were obtained using a simulation domain for the fluid model of 6 km length and 2 km radius. In this case, 2 × 10 15 electrons are injected into an electric field corresponding to 3.57 kV/cm (or 1.27 RREA thresholds) at ground level, but scaled to 12 km altitude. The validity of these assumptions will be discussed in Comparison of a Gaussian fit of our simulated current propagated using an FDTD model (black), and the analytical method outlined in Section 2.2.1 (orange and purple). The orange line represents a case with no ground attenuation applied (conducting ground), while the purple line corresponds to ground attenuation based on propagation over land. The FDTD model is using ground attenuation based on propagation over seawater, which, being more conductive than land, makes the result fall between the two analytical cases.

Section 4. The initial field extends over the first two thirds of the 6 km domain to avoid interference from the simulation boundary. The background field is initially set as zero in the last third, so as to ensure the dampening of the RREAs before they reach the end of the simulation domain. The electrons are injected following a Gaussian distribution in time with the duration of the injection being t inj = 15 μs, corresponding to a standard deviation of 10.6 μs.

For the sake of illustration, Figure 6 shows the electron density obtained in the model at t ≃ 50 μs, for which moment of time it reaches a maximum of ∼10 14 m -3 .

Figure 7 shows a map of the total current obtained using these initial conditions. The total current is obtained through integration of the current flux over cross-sections of the simulation domain. The total current flux is the sum of positive and negative ions, electrons, and the displacement current. The z-axis shows the current along the length of the simulation domain. For a resolution Δz ≃ 1 m, Figure 7 can be seen as a map of current moment elements, which relates it directly to Equation 14. The abscissa shows the development in time. The duration of injection of 15 μs determines the rise time of the current, which peaks at 74.2 kA. The decay of the current is primarily driven by three-body attachment. The remaining tail of current is due to ion drift in the electric field, and falls off slowly over a timescale of several milliseconds. In 3 ms, the tail reduces to 60% of its amplitude right after the initial peak (i.e., ∼30 μs after the initial peak).

Figure 8 shows the resulting radio signal as it would be registered by the detectors used by [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF]. The signal has been propagated over a ground distance of 156 km from the source using the analytical formulation described in Section 2.2.1. Filters have then been applied to approximate the filter response (Section 2.3). The result is compared with data collected at the FT, on 4 September 2015 [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF].

As can be seen, the simulation result and the slow LF pulse are similar in shape, amplitude, and duration. The asymmetry of the first and second positive peaks is present in both cases, with a ratio of the first to second positive peak being 1.63 for the slow LF pulse and 1.70 in the simulation. The difference in amplitude between the positive and negative peaks is a bit larger-the ratio between the first positive peak and the negative peak is -0.74 for the slow LF pulse, and -0.58 for the simulation. The duration of the slow LF pulse is given by [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF] as 65 μs. The duration of the current obtained by our simulation is similar; ± 2 ⋅ t inj is 60 μs. Our simulation results have a slightly slower decrease toward the end of the signal than the slow LF pulse. The injection time of the seed electrons is a free parameter in our simulations, and has as such been chosen to fit the slow LF pulse as closely as possible.

Discussion

Validity of Initial Conditions

The context of the simulations performed in this paper is that of a lightning leader injecting seed runaway electrons into a large-scale, homogeneous electric field. The number of seed electrons used varies around 10 15 , which is on the lower side, but consistent with previous papers using such a context (e.g., [START_REF] Celestin | Variability in fluence and spectrum of high-energy photon bursts produced by lightning leaders[END_REF]Mailyan et al., 2019). The electric fields we use are, with small variations, around 1.3 times the threshold RREA field [START_REF] Dwyer | High-energy atmospheric physics: Terrestrial gamma-ray flashes and related phenomena[END_REF]. Since no RREAs will happen below the threshold field, this is then a relatively low field to sustain RREA-related processes. We assume the air density of 12 km altitude, and this is also where we place the peak of the current when propagating its radio signal. This is consistent with previous literature on TGFs (e.g., [START_REF] Cummer | The source altitude, electric current, and intrinsic brightness of terrestrial gamma ray flashes[END_REF]Xu et al., 2012), and it is also within the range of source altitudes for the slow LF pulses as described in [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF]. Finally, using these initial conditions, we find the number of runaway electrons at the midpoint of our simulation domain to be in the 10 17 to 10 18 range (see Table 2), which is consistent with the number of runaway electrons needed to explain TGF fluences at satellite altitude (e.g., [START_REF] Celestin | Variability in fluence and spectrum of high-energy photon bursts produced by lightning leaders[END_REF].

Comparison to Other Observations

In Figure 9, the results presented in Section 3 are compared to two other detections of slow LF pulses from FT. While the first observation was done 156 km from the source, these two are significantly farther away. One was recorded on 3 August 2010, at a distance of 466 km from the source, and the other on 10 August 2014, 472 km from the source. Our simulation result is clearly a poorer match for these events, noticeably the amplitude of the signal is significantly lower than what was measured. This implies that the different observations presented in [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF] are inherently different-these two events correspond to a stronger source signal than the event seen in the previous section. We conclude that there must be a range of initial physical conditions capable of leading to slow LF pulses. 8, here propagated to 466 km (left) and 472 km (right) to match detections done at FT on 3 August 2010 and 10 August 2014, respectively. In the observational data (red) one can see the ionospheric reflection of the signal following the slow LF pulse, which is not included when using our analytical approach to radio propagation. BERGE ET AL. Examples of better fits for these two events can be seen in Figure 10. Case (a) is the observation from 2010, which was 466 km away from the source. Case (b) is from 2014, at 472 km distance. Due to the significant computation time, no automatic fitting procedure is used here and the initial parameters of the simulations for these cases are summarized in Table 1. The simulation presented in Section 3 is also included for comparison. The simulation domain length and radius remain 6 and 2 km, respectively, through all the cases. The external electric field extends for the first 2/3 of the length throughout. The peak current and the number of runaway electrons at the mid-point of the simulation are included in Table 2. The peaks following the slow pulses here are due to ionospheric reflection of the signal. These results are not necessarily unique, that is, there might be other configurations of initial parameters that give equally good matches. The following section will describe the impact of varying the free parameters in our model.

Contributions of Different Initial Parameters

Different initial conditions of the simulation will naturally yield different results. The free parameters we have in the current model include the duration of injection, the number of initial seed electrons, the amplitude of the external electric field, the extent of the electric field, and the radius of injection. The time-profile of electron injection follows a Gaussian function throughout the paper. Interestingly, we find that moderate changes of the initial RREA radius produce negligible changes in the final result. The injection time has been chosen to fit the observations we have access to. This parameter controls both the shape and amplitude of the result, as well as of course its duration. Changing the initial number of electrons and the initial electric field will also affect the results' shape and amplitude. In general, more electrons, a stronger field, and a larger extent of the field will all lead to a stronger resulting pulse. This also means that decreasing some and increasing others of these values to some degree can balance each other out, creating a range of values that can give reasonable results. Since the parameters all affect each other, the number of possible combinations quickly becomes large, and we have not performed a full analysis of the whole parameter space. Our model also has numerical limitations. In particular, when the initial conditions lead to very high particle densities, the time step of the simulation becomes too small to get results within any reasonable time frame.

To give the reader a feel for how changing the initial parameters may change the results, we present some examples in Figure 11. We have used the case presented in Section 3 as a reference, and varied one parameter at a time. As can be seen, even small changes in the initial parameters produce visible changes of both the shape and amplitude of the signal. Increasing the injection time increases the duration of the signal as expected, and also diminishes the amplitude. Increasing the electric field, number of seed electrons, and the length of the simulation domain (and thereby the extent of the electric field, which is still two-third of the domain length) all cause an increase of the amplitude of the signal, albeit in subtly different ways. The shape of the signal, that is, the ratios between the different peaks, changes in different ways for the different parameters. The peak current and resulting number of runaway electrons are also different. A summary of these properties of the results can be seen in Table 2.

The variations shown in Table 2 and Figure 11 are small, and cannot be considered physically less reasonable choices for initial parameters than those of the reference case from Section 3. Combined with the inherent difference of the three observations of slow LF pulses presented in Sections 3 and 4.2, this suggests that a range of physical conditions can and do produce LF pulses with recognizable characteristics. While all the parameters presented cause a change in the final signal, future refinement of the simulations and larger data sets of slow LF pulses may allow for further constraints on the physical conditions necessary for slow LF pulse production. The table value is the electric field scaled to ground-level. The simulation is performed at 12 km altitude, which is consistent with a source altitude of 10-15 km found by Pu et al. (2019). b Corresponds to the number of runaway electrons in the mid-plane of the simulation domain. 11. The first column contains the values of the reference case (the solid black curve in Figure 11), while the four following columns correspond to the variations presented in dashed grey lines in the four subplots of Figure 11.

Conclusions

1. Using a multifluid plasma model assuming a number and time scale of seed electrons as an initial condition, we find that TGF-producing RREAs in the framework of leader-based models are consistent with the emission of slow LF pulses 2. The timescale of the final signal is dependent on the timescale of the initial electron injection. Using an injection time of 15 μs, we get a final signal of similar duration as observations of slow LF pulses [START_REF] Pu | Low frequency radio pulses produced by terrestrial gamma-ray flashes[END_REF] 3. The magnitude of the electric field, the length of the simulation domain, the injection time of seed electrons, and the initial number of electrons all change the shape and amplitude of the signal. This constrains the combinations of parameters which yield results in agreement with observational data 4. This study opens the question as to what underlying processes lead to the ∼15 μs injection timescale of runaway electrons The black line in all four subplots is the case presented in Section 3. The gray dashed lines each represent a variation of this case, where the value of one parameter has been changed slightly, while all the others have been kept the same as the case in black. That is, in the first subplot only the electric field has been changed, not the parameters presented in the other subplots.

Figure 1 .

 1 Figure 1. Schematic representation of the propagation of information in the model. Seed electrons are injected according to an injection function, I(t) (red). For each time step, δt, corresponding to the time it takes the RREA to move across one simulation grid size, Δz, new values for the diffusion radius, R D , and exponential term, ξ, are calculated. The variables are dependent on their previous values, and as such at the current time of the simulation (shaded light red), each grid cell uses information from where the corresponding electrons were at previous iterations (following the diagonal lines).

Figure 2 .

 2 Figure 2. Sketch of the geometry of the antenna above a perfectly conducting ground, whose emission is described by Equation 14. Based on Figure 2 in Uman et al. (1975).

Figure 3 .

 3 Figure3. A diagram representing the filters applied to our simulation results to approximate ground attenuation (first group of filters, left), and to mimic the detector response of the LF sensor at FT (second group of filters, right). The small gray boxes represent scaling factors applied to the signal. All Butterworth filters are first order, while the Chebyshev filters have the order noted. The single pole 1-MHz Butterworth filter acts as a differentiator, which is part of the detector response. For ground attenuation at distances not mentioned here, the distance is rounded to the closest of these three values for the purpose of ground attenuation application only.

Figure 4 .

 4 Figure 4.Comparison of a Gaussian fit of our simulated current propagated using an FDTD model (black), and the analytical method outlined in Section 2.2.1 (orange and purple). The orange line represents a case with no ground attenuation applied (conducting ground), while the purple line corresponds to ground attenuation based on propagation over land. The FDTD model is using ground attenuation based on propagation over seawater, which, being more conductive than land, makes the result fall between the two analytical cases.

Figure 5 .

 5 Figure 5. Schematic representation of the information flow during modeling. The black rectangles represent steps of the modeling process. The input into the next step is the output from the previous, shown in gray squares. This figure provides a general overview, for details concerning the ground attenuation and detector response, see Figure 3.

Figure 6 .

 6 Figure 6. Electron density obtained in the model at t ≃ 51 μs. The initial number of runaway electrons is 2 × 10 15 under an electric field of 1.27 RREA thresholds. The upper panel presents a cross-sectional view of the electron density while the lower panel shows the electron density along the z-axis at the same time in a semi-logarithmic plot.

Figure 7 .

 7 Figure 7. Map of the current produced by our fluid simulation along the z-axis, and its evolution in time.

Figure 8 .

 8 Figure 8. Simulation results compared to slow LF pulse. The LF pulse was detected 156 km from the source, and the simulated radio emission has been propagated the same distance. The curves have been aligned in time.

Figure 9 .

 9 Figure9. The same simulation as shown in Figure8, here propagated to 466 km (left) and 472 km (right) to match detections done at FT on 3 August 2010 and 10 August 2014, respectively. In the observational data (red) one can see the ionospheric reflection of the signal following the slow LF pulse, which is not included when using our analytical approach to radio propagation.

Figure 10 .

 10 Figure 10. Simulations of the current created by the RREA (left), and the resulting radio pulse (right). The initial parameters of the simulation are made so that the radio pulse (black, right) fits the observations (red, right) from 3 August 2010 (a) and 10 August 2014 (b). See text for description of the parameters. In the observational data, one can see the ionospheric reflection of the signal following the slow LF pulse.

Figure 11 .

 11 Figure 11. Variation of different initial parameters of the simulation. The black line in all four subplots is the case presented in Section 3. The gray dashed lines each represent a variation of this case, where the value of one parameter has been changed slightly, while all the others have been kept the same as the case in black. That is, in the first subplot only the electric field has been changed, not the parameters presented in the other subplots.

Table 1

 1 Table of Initial Parameters, as Well as Some Selected Results, for the Three Events Note. These are numerical values for the results presented in Figure

		Reference	E = 3.60 kV/cm	N = 3 ⋅ 10 15	t inj = 17 μs	l = 6.1 km
	Peak current (kA)	74.24	93.07	93.38	66.68	86.63
	Number of runaway electrons	7.10 × 10 17	9.00 × 10 17	9.01 × 10 17	7.12 × 10 17	8.35 × 10 17
	Ratio of positive peaks	1.70	1.66	1.84	1.71	1.60
	Ratio of first peak to negative peak	-0.58	-0.59	-0.60	-0.58	-0.58

Table 2

 2 Table of Selected Results While Varying the Initial Parameters One by One, Compared to the Reference Case Presented in Section 3

Acknowledgments

Nini L. Berge, Matthieu B. Garnung, and Sebastien Celestin's research was funded by the French Space Agency (CNES) through projects OREO and TARANIS, and by the French Region Centre-Val-de-Loire. Sebastien Celestin also acknowledges support from Institut Universitaire de France (IUF). Wei Xu and Robert A. Marshall acknowledge support from the National Science Foundation (NSF) grant AGS-1732359. Steve A. Cummer acknowledges support from the NSF through grant AGS-2026304.

Data Availability Statement

The simulation results used to produce the figures and the associated analysis are available at https://doi. org/10.6084/m9.figshare.c.5660416.v1.