Determining Atmospheric Electric Fields using MGMR3D

T. N. G. Trinh [email protected] Physics Education Department, School of Education, Can Tho University, Campus II, 3/2 Street, Ninh Kieu District, Can Tho City, Viet Nam O. Scholten [email protected] Kapteyn Astronomical Institute, University of Groningen, P.O. Box 72, 9700 AB Groningen, Netherlands Interuniversity Institute for High-Energy, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium S. Buitink Astrophysical Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium Department of Astrophysics/IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands K. D. de Vries Astrophysical Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium P. Mitra Astrophysical Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium Faculty of Physics, University of Warsaw - Warsaw, Poland T. Phong Nguyen Department of Physics, Can Tho University, 3/2 Street, Ninh Kieu, Can Tho City 94000, Vietnam D. T. Si Department of Education and Training, 3/2 Street, Ninh Kieu District, Can Tho City, Viet Nam

Abstract

Cosmic-ray particles impinging on the atmosphere induce high-energy particle cascades in air, an Extensive Air Shower (EAS), emitting coherent radio emission. This emission is affected by the presence of strong electric fields during thunderstorm conditions. To reconstruct the atmospheric electric field from the measured radio footprint of the EAS we use an analytic model for the calculation of the radio emission, MGMR3D. In this work we make an extensive comparison between the results of a microscopic model for radio emission, CoREAS, to obtain an improved parametrization for MGMR3D in the presence of atmospheric electric fields, as well as confidence intervals. The approach to extract the electric field structure is applied successfully to an event with a complicated radio footprint measured by LOFAR during thunderstorm conditions. This shows that, with the improved parametrization, MGMR3D can be used to extract the structure of the atmospheric electric field.

cosmic rays; shower parameters; atmospheric electric fields; radio emission; extensive air showers

I Introduction

Lightning is a common phenomena but the detailed understanding of its generation and development is still unknown. One of the reasons why this topic is still under investigation is that atmospheric electric fields inside thunderclouds are difficult to measure. Existing measurements, from balloons and airplanes, are limited because the measurements depend on the paths of the balloons or the aircrafts will affect the nature of the thunderclouds. A new method to determine atmospheric electric fields is using their influence on radio emission emitted from extensive air showers [1]. Unlike the balloon and airplane measurements, this unique tool is not limited by the wind conditions and is sensitive to a large part of the cloud.

When a cosmic ray enters the atmosphere of the earth, it interacts with air molecules and generates a particle cascade, called an extensive air shower (EAS). The particles in the EAS move with velocities near the speed of light and are concentrated in the thin shower front. In an EAS during fair-weather conditions, called fair-weather showers, the electrons and positrons are deflected in opposite directions by the Lorentz force caused by the geomagnetic field. This induces an electric current in the shower front that is transverse to the shower axis. Since the electric current changes as a function of height, due to the change of the number of particles in the EAS, radio waves are emitted [2, 3, 4]. A secondary, yet important, contribution to the radio emission is caused by an excess of electrons in the shower front due to knock-out from atmospheric molecules by shower particles. This process creates a radio pulse that is linearly polarized but oriented radially to the shower axis [5, 6]. This charge excess emission interferes with the geomagnetic emission that is linearly polarized along the direction of the Lorentz force.

An EAS that occurs during thunderstorm conditions when there are strong atmospheric electric fields, is called a thunderstorm shower, and produces radio emission that differs considerably from that of fair weather showers [1, 7, 8]. In thunderstorm showers, the atmospheric electric field exerts an electric force on the electrons and positron that is usually much stronger than the Lorentz force and this affects the radio emission. The component of the electric field perpendicular to the shower axis will change the direction and magnitude of the transverse current changing the shape of the radio footprint. The electric field parallel to the shower axis accelerates the electrons or positrons, depending on its direction, and thus they gain more energy. As a result, these particles will generate additional low-energy particles. However, these additional particles travel with a velocity smaller than the speed of light and thus they trail far behind the shower front. Thus, their radiation is not added coherently in the frequency range from 30 MHz to 80 MHz at which typical cosmic-ray airshower radio detectors operate [7].

There are several models to describe radio emission from air showers. Microscopic models such as ZHAires [9] and CoREAS [10] are based on full Monte-Carlo simulation codes. Macroscopic models such as MGMR [2], EVA [3] and Radio Morphing [11] calculate the emission of the bulk of electrons and positrons described as currents. SELFAS2 [12] is a mix of macroscopic and microscopic approaches. Recently, MGMR3D [13] has been introduced that uses a multi-dimensional parametrization of the current density in the air showers to calculate the radiation field using Maxwell equations.

Because atmospheric electric fields influence the radio emission from air showers, the radio footprint, as is measured on the ground, can be used to determine the strength, direction, and the altitude dependence of the atmospheric electric field along the path of the air showers [14]. Since we have only a forward model, i.e. calculate the footprint by assuming an electric field configuration, we have to apply a search algorithm to extract the field structure. This only works with a fast model calculation that is deterministic, i.e. non Monte-Carlo based. MGMR3D is potentially such a model. However, MGMR3D relies on a parameterization of the emitting currents. To gauge these parameters we perform an extensive comparison between MGMR3D and the full Monte-Carlo code CoREAS.

In Ref. [13, 8] a parametrization was proposed based on CONEX-MC [15] simulations. In particular it had been assumed that the induced drift velocity, driving the currents, is inversely proportional to the air density, $v_{d}\propto 1/\rho$ . This assumes that the drift velocity is inversely proportional to the friction force, the collision frequency. In Ref. [16] it was shown that an inverse dependence of the drift velocity on $v_{d}\propto 1/\sqrt{\rho}$ yields a much better agreement with the results of microscopic calculations. The $1/\sqrt{\rho}$ dependence is characteristic of an equilibrium situation, terminal velocity, where the acceleration of the fast electron by the Lorentz force is compensated by the average friction due to collisions with air molecules [16]. This has led to a considerable improved parametrization for the description of fair weather showers in MGMR3D [16]. This different proportionality for the drift velocity has made us to revisit the parametrizations that enter in MGMR3D in the presence of atmospheric electric fields, by making an extensive comparison with a large number of CoREAS calculations for random shower directions and various structures of the atmospheric electric fields. The improved parametrization allows us to extract the fields for more complicated configurations that were inaccessible before.

In Section II we present a very short review of the most important aspects of CoREAS that describes radio emission on a microscopic level and which we consider as ‘the truth’ to which we gauge the parameters that enter in the modeling of the currents entering in MGMR3D, as discussed in Section III. We first review the essentials of the parametrization for fair-weather showers [16] after which we present the discussion of the parametrization in the presence of electric fields. In Section IV we apply the improved model to data for the event which was not possible to analyze in Ref. [8] to extract the atmospheric electric field.

II CoREAS simulations

To gauge the parameters in the MGMR3D model, discussed in Section III, we perform extensive comparisons with the radio emission footprints from EAS as calculated from CoREAS, which we regard as the Monte-Carlo truth. CoREAS is a plug-in for the particle simulation code CORSIKA [17] and is based on a microscopic description of the radiation mechanism. In CoREAS the radio emission from each electron and positron in the shower is calculated, without making any assumptions on the type of radiation. In other words, CoREAS does not simulate the transverse-current and charge-excess components separately, but it produces the complete radiation field from the emissions of the individual particles. Therefore, CoREAS simulations are compute intensive with running times of the order of days. Using the electric field EFIELD option [18] in CORSIKA it is possible to simulate air showers passing through strong electric-field regions. For this work the radiation profile is calculated in the shower plane (the plane perpendicular to the shower axis passing through the point where the shower axis touches ground), with the $\hat{x}$ -axis in the direction of ${\bf{v}\times\bf{B}}$ , along the direction of the Lorentz force, and $\hat{y}$ -axis along ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ . Here $\mathbb{v}$ is the velocity of the shower and $\mathbb{B}$ is the geomagnetic field. We use the shower plane instead of the ground plane since it allows for an easier interpretation. To construct a two-dimensional radio footprint, we run simulations for 160 antenna positions in a star-shaped pattern in the shower plane with eight arms, where each arm contains 20 antennas with a spacing of 25 m [19].

To compare with the results from MGMR3D, we extract the values for the Stokes parameters at the 160 antenna positions and thus capture both intensity and polarization of the radio emission. Because the aim of this work is to develop a realistic approach to extract electric fields from LOFAR [20] data, we construct the Stokes parameters from the CoREAS simulations in the frequency band from 30 MHz to 80 MHz. For each position in each simulated shower, the Stokes parameters are calculated as

I=\frac{1}{n}\sum_{i=0}^{n-1}\left(\left|\varepsilon_{i,{\bf{v}\times\bf{B}}}\right|^{2}+\left|\varepsilon_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right|^{2}\right),

(1)

Q=\frac{1}{n}\sum_{i=0}^{n-1}\left(\left|\varepsilon_{i,{\bf{v}\times\bf{B}}}\right|^{2}-\left|\varepsilon_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right|^{2}\right),

(2)

U+iV=\frac{2}{n}\sum_{i=0}^{n-1}\left(\varepsilon_{i,{\bf{v}\times\bf{B}}}\varepsilon^{*}_{i,\bf{v}\times\left({\bf{v}\times\bf{B}}\right)}\right).

(3)

$\varepsilon_{i,x}$ is the complex-valued signal radiation field along direction $x$ , where $i$ denotes the sample number and $x={\bf{v}\times\bf{B}}$ or $x=\bf{v}\times\left({\bf{v}\times\bf{B}}\right)$ . $n$ is the number of time samples centered around the pulse-peak position. Stokes $I$ is the intensity of the radio emission. Stokes $Q$ and $U$ are used to derive the linear polarization angle

\psi=\frac{1}{2}\mathrm{tan}^{-1}\left(\frac{U}{Q}\right),

(4)

and Stokes $V$ is the intensity of the circularly polarized fraction of the radiation.

In order to perform chi-square fitting with realistic sensitivities, error bars (confidence levels) were assigned to the values of the Stokes parameters for the CoREAS simulated showers, given by

\sigma_{k,\left(I,Q,U,\mathrm{or}\,V\right)}=2\sqrt{\frac{\xi_{I}}{n}(I_{k}+\xi_{I})},

(5)

where $n$ is the number of time samples as enters in the equations for the Stokes parameters Eq. 3, $I_{k}$ is the value of Stokes $I$ for a test antenna at position $k$ and $\xi_{I}$ is taken at the level of the instrumental noise measured in LOFAR antennas [19]. This definition of the error bars of the Stokes parameter guarantees that the relative error bars at positions near the core, where the intensity is usually large, is small while at large distances where the intensity is small, the relative error bars become large.

III Parameterization in MGMR3D

In MGMR3D [13], the Maxwell’s equations are solved in the far-field approximation to obtain the radio footprint, using a parametrized charge-current density that depends on atmospheric depth, radial distance from the shower axis, and distance from the shower front. In earlier work, presented in Ref. [8], the parametrization of the charge-current density was obtained through a comparison of selected cases with the results of CONEX-MC and CoREAS calculations. The resulting parametrization was adequate for extracting the structure of relatively simple atmospheric electric fields, however for a thunderstorm EAS with a complicated radio footprints there were are large discrepancies between the MGMR3D and CoREAS results [8].

Following the approach presented in Ref. [16] we tune the parametrization in MGMR3D by making a statistical comparison of the calculated radio footprint with the results of microscopic calculations for a large numbers of showers passing through different electric-field structures. The advantage of such an approach is that a comparison is made at the level of the important observables for later applications. In addition this approach allows to asses the intrinsic accuracy of obtaining the structure of the atmospheric electric field.

Since an extensive discussion of the parametrization for fair-weather showers is presented in Ref. [16], we will in this paper only summarize the essentials of the parametrization for fair-weather showers in Section III.1. The results [16] show that a good agreement between MGMR3D and CoREAS is obtained for fair-weather showers. The parametrization of the charge-current cloud in the presence of electric fields is discussed in Section III.2.

III.1 Parameterization for fair-weather showers

Under fair-weather conditions only the Lorentz force, $\vec{F}_{\perp}=e\vec{v}\times\vec{B}$ , is acting on the charged particles in the EAS, causing the electrons and positrons to drift in opposite directions and thus induce an electric current. As presented in [16, 21] the induced drift velocity depends on the square-root of the air density and can be parametrized as

\vec{v}_{\perp}(h)=\vec{F}_{\perp}/F_{0}\frac{\sqrt{\rho(X_{\rm max})/\rho(h)}}{((X_{\rm max}-X_{v})/(X_{z}(h)-X_{v})+3.)}\;,

(6)

with $X_{v}=50$ g/cm² and $F_{0}=62.5$ keV/m. This -on first sight- surprising dependence on $\sqrt{\rho}$ is due to the fact that at the equilibrium sideways ‘friction’ force, proportional to $\rho v_{\perp}^{2}$ , equals the transverse component of the Lorentz force. Since only the Lorentz force is acting, the drift velocity $\vec{u}_{\perp}$ is still non-relativistic and it is parameterized as

\vec{u}_{\perp}(h)=c\vec{v}_{\perp}(h)/\sqrt{1+v_{\perp}^{2}/u_{0}^{2}}\;,

(7)

where the parameter $u_{0}$ is taken much larger than typical fair-weather drift velocities, resulting in a linear scaling of the drift velocity with the Lorentz force.

In MGMR3D the number of charged particles at a radial position $r_{s}$ from the shower axis and a longitudinal position $d_{h}$ from the shower front is written as $Nf(d_{h})w/r_{s}$ , with $w(r_{s})$ the normalized radial particle distribution, and $f(d_{h})$ the normalized longitudinal particle distribution. The current is the product of the charged particle number and the drift velocity, where both depend on the penetration depth in the atmosphere, $X_{z}(h)$ . The current is thus expressed as

\vec{J}_{\perp}(t_{s},x_{s},y_{s},h)={w(r_{s})\over r_{s}}\,f(d_{h},r_{s})\,N_{c}(X_{z})\vec{u}_{\perp}(h)\;,

(8)

at a radial distance $r_{s}=\sqrt{x_{s}^{2}+y_{s}^{2}}$ from the shower axis. The parametrization of the longitudinal shower profile for the current in the shower front is based on the Gaisser-Hillas formula [22] for the dependence of the number of charged particles on $X_{z}$ ,

N_{c}(X_{\rm z})=\left({X_{\rm z}-X_{0}\over X_{\rm max}-X_{0}}\right)^{X_{\rm max}-X_{0}\over\gamma}e^{X_{\rm max}-X_{\rm z}\over\gamma}

(9)

where $\gamma$ is a parameter controlling the width of the distribution and $X_{0}$ the first-interaction point. $X_{\rm max}$ is the atmospheric depth where the number of secondary particles reaches a maximum in the air shower.

The spatial distribution of the charged particle density is parametrized as

w(r_{s})=N_{w}\,\xi(\xi+1)^{-2.5}\;;f(d_{h},r_{s})=N_{f}\,{\eta\over e^{\sqrt{\eta}}+1}\;,

(10)

with $\eta={d_{h}/\lambda}$ where $d_{h}$ is the distance from the shower front, and $\xi=r_{s}/R_{0}$ is a scaled radial distance where we introduced the radiation radius $R_{0}$ . The scaling factors $N_{w}$ and $N_{f}$ are chosen such that the distributions are normalized, $\int w(r)\,dr=1$ and $\int_{0}^{\infty}f(d_{h},r_{s})dd_{h}=1$ .

It was observed that the optimum value for $R_{0}$ depends on the distance from $X_{\rm max}$ to the shower impact point. For small distances there is an almost linear dependence until a saturation value of $R_{0}=50$ m is reached at a distance of about 5 km, independent of zenith angle. A full account of the parametrization for fair weather showers in MGMR3D is presented in [16].

III.2 Parameterization for thunderstorm showers

Due to convection flows in clouds there is a build up of electric charge that has a layered structure when considering large clouds. In some cases this can spontaneously cause a lightning discharge. The charge layer structure depends on detailed cloud dynamics and temperatures, see for a recent reference [23]. To extract the atmospheric field configuration we thus parameterize a semi-realistic atmospheric field configuration as built in three different layers where in each layer the electric field may have a different orientation and strength. Including the boundary heights we thus arrive at nine parameters for the field in addition to two parameters for the location of the shower core as well as one for the energy of the cosmic ray and one for $X_{\rm max}$ . We assume that the arrival direction of the cosmic ray is known from the arrival time of the radio signals in the different antennas. Since this adds up to a large number of parameters it is imperative to limit the introduction of additional parameters that are not likely to greatly influence the structure of the extracted electric fields. The values for the radiation radius $R_{0}$ , the amount of charge excess, parameters defining the structure of the shower such as the first interaction depth $X_{0}$ , and the width of the shower profile $\lambda$ are therefore kept fixed to the generic values used for fair-weather showers. The values of these parameters are given in Table. 5 in the Appendix (except for $R_{0}$ that was discussed in the previous section).

There are two parameters that are more specific for the cases in which there is an electric field present. One is related to how quickly the induced transverse currents adjust to an electric field that changes with height, discussed in Section III.2.1. Another parameter is $u_{0}$ as introduced in Eq. (7). In Section III.2.2 the optimal values for these additional parameters are obtained by fitting an ensemble of showers simulated in CoREAS, following the approach introduced in [16].

To verify the implementation of atmospheric electric fields in MGMR3D and to test the superposition principle we present in Section III.2.3 a comparison with CoREAS calculations is performed for more intricate field configurations.

III.2.1 Height-dependent electric fields

In the presence of atmospheric electric fields the net force acting on the particles changes. The component of the electric field parallel to the shower axis is assumed to be well below the runaway breakdown [24]. This component can increase the number of secondary non-ultra-relativistic (with energies well below a GeV) electrons while decreasing the number positrons or the other way around, depending on the its orientation. However, these particles are generally non relativistic and thus they trail far behind the shower front [7]. For this reason, they do not contribute much to the radio radiation in the frequency range from 30 MHz to 80 MHz [7] and thus this field component can be ignored.

We have noticed that electric fields in excess of 50 kV/m, that are transverse to the shower, do affect the structure of the shower profile. We see an increase in the number of charged particles starting at the height where the electric field perpendicular to the shower axis, $\vec{E}_{\perp}$ , is present. As a result, the value of $X_{\rm max}$ is modified in the presence of a strong electric field. The reason for this effect is still unknown. To model this we modify the parametrization of the shower profile as

\displaystyle N_{c,E}(X_{\rm z})=N_{c}(X_{\rm z})\left(1+0.0015\,\frac{\rho(0)}{\rho(h)}\frac{F_{\perp}^{2}}{|v\times B|^{2}}\right)\;,

(11)

where $F_{\perp}$ is given in Eq. (12).

As has already been argued, $\vec{E}_{\perp}$ , has a large effect on the radio emission. In the presence of atmospheric electric fields the transverse force acting on the shower particles is

\vec{F}=q(\vec{E}_{\perp}+\bf{v}\times\bf{B})\;.

(12)

This changes both the magnitude and the direction of the drift velocity which is parallel to $\vec{F}$ . In strong atmospheric electric fields, the net forces become large and thus the drift velocity increases. The simulations show that the radio emission does not increase proportionally and shows saturation effect for strong fields [7]. This can be understood from relativistic effects where the energy is no longer proportional to the square of the velocity. To account for this we follow Ref. [7] and introduce the saturation effects in the drift velocity, Eq. (7), by means of the parameter $u_{0}$ . In Section III.2.2, Table 2, the value for this parameter is discussed.

Atmospheric electric fields may vary strongly with height. Thus another important parameter is the time, or distance, it takes for the drift velocity to adapt to the changes in the electric field. To parametrize this we replace $\vec{F}_{\perp}$ in Eq. (6) by an effective force, $\vec{\tilde{F}}_{\perp}$ , defined as

	$\displaystyle\vec{\tilde{F}}(h)$	$\displaystyle=$	$\displaystyle\vec{F}(h)+$		(13)
			$\displaystyle\sum_{j}\frac{1}{1+e^{(a+(h-H_{s}-h_{j})\,b)}}e\Delta\vec{E}_{j}\;,$		(13)

where

a=D_{Es}\,H_{s}\,\rho_{\rm ground}/X_{EM}\;,

(14)

and

b=D_{Es}\,\rho(h_{j})/X_{EM}\;.

(15)

Here $\Delta\vec{E}_{j}$ is the change in the electric field at height $h_{j}$ and $X_{EM}$ is the interaction length for an electromagnetic particle. The parameter $D_{Es}$ (for E-smooth) governs the distance over which the drift velocity (through the effective force) adapts to the change in the electric field, if $D_{Es}$ is large the second form will reduce to a step function $\theta(h+H_{s}-h_{j})$ . The drift velocity will only adapt to the applied force after a certain relaxation time and the parameter $H_{s}=700$ is chosen such that the average shift in height, as shown later in Fig. 2 is minimal.

Refer to caption — Figure 1: The ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ -component of the current profile for $D_{Es}=4,6,8$ and $10$ as a function of atmospheric depth.

Fig. 1 shows how the current along the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ direction changes when the EAS with $X_{\rm max}=500$ g/cm² passes through an electric field boundary for different values of $D_{Es}$ . An electric field of 50 kV/m is taken between a height of 4.5 km ( $X=590$ g/cm²) and the ground, oriented along the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ -direction. It shows that the current changes faster when the value of $D_{Es}$ is larger. The drop-off of the current at large values of $X$ , i.e. close to the ground, is due to the decrease in the number of shower particles.

To determine the optimal values for the parameters $D_{Es}$ and $u_{0}$ we simulate a set of 94 vertical showers passing through a one-layer electric field using the microscopic code CoREAS. For each shower there is an electric field present from a certain height $h$ to the ground with certain strength $E$ and angle $\alpha$ where $\alpha$ is the angle between the electric field and ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ -direction. The height $h$ , the strength $E$ , the angle $\alpha$ , as well as $X_{\rm max}$ are chosen randomly.

To select the optimal value for $D_{Es}$ , we choose four different values for $D_{Es}$ , 4, 6, 8, and 10. As shown in Table. 1, the mean value of $\chi^{2}$ does not change much when changing $D_{Es}$ and thus we do not need to determine $D_{Es}$ to higher precision. While keeping this parameter fixed and using the height $h$ , the strength $E$ and the angle $\alpha$ obtained from CoREAS, we perform a chi-square fit of $X_{\rm max}$ . We also keep $H_{s}$ fixed at 700 and $u_{0}=10$ . By using a Levenberg-Marquardt minimization procedure based on a steepest descent method, we optimize the value of $X_{\rm max}$ by minimizing

\chi^{2}=\sum_{\text{antenna}}\sum_{S=I}^{Q,U,V}\left(\frac{S_{k,C}-f_{r}S_{k,3D}}{\sigma_{k}^{2}}\right)\,,

(16)

where $S_{k,C}$ are the Stokes parameters from CoREAS calculations for antenna at position $k$ and $S_{k,3D}$ are the ones from MGMR3D. $\sigma_{k}$ is the uncertainty as defined in Eq. (5) and $f_{r}$ is the normalization factor for the radio intensity. The results indicate that the value of $D_{Es}$ has only a minor effect on the fitted values of $X_{\rm max}$ . However, as shown in Table. 1, the $\chi^{2}$ has a shallow minimum for $D_{Es}=8$ . Therefore, we keep $D_{Es}$ fixed to this value for all subsequent calculations. The results of fitting also show that the value of $\chi^{2}$ does not depend on $X_{\rm max}$ and the electric field configuration.

$D_{Es}$	4	6	8	10
Mean $\chi^{2}$	0.183	0.178	0.175	0.178
Standard deviation	0.156	0.150	0.146	0.145

Table 1: Mean

\chi 2

and its standard deviation for different values of

DEs

The other parameter which is kept fixed for all thunderstorm showers is $u_{0}$ as introduced in Eq. (7). Following the same procedure as for $D_{Es}$ , fitting $X_{\rm max}$ while keeping $u_{0}$ , $h,E$ and $\alpha$ fixed, the mean $\chi^{2}$ and the standard deviation are calculated for several values of $u_{0}$ in the range from 0.3 to 10. We also see that the parameter $u_{0}$ hardly influences the extracted values of $X_{\rm max}$ . As shown in Table 2, the mean value of $\chi^{2}$ decreases when $u_{0}$ increases. However, when $u_{0}$ is equal or larger than 1, it does not affect the quality of the fit anymore, so we keep it fixed at the value of 10.

$u_{0}$	0.3	0.6	1	3	7	10
Mean $\chi^{2}$	0.188	0.177	0.176	0.176	0.175	0.176
Standard deviation	0.164	0.149	0.147	0.146	0.147	0.146

Table 2: Mean

\chi 2

and its standard deviation for different values of

u_{0}

III.2.2 Accuracy of determined electric fields

We first investigate the accuracy of extracting one-layer electric fields. For all 94 vertical showers passing through a one-layer electric field simulated by CoREAS as discussed in the previous section, we follow a fitting procedure using MGMR3D to reconstruct the structure of the field and compare this with the true values used in the CoREAS calculation. In this section the true electric field is taken homogeneous from a top height $h_{\rm true}$ to the ground with a strength $E_{\rm true}$ , making an angle $\alpha_{\rm true}$ with the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ axis. The reconstructed values, obtained by fitting the radio footprint using MGMR3D, carries a subscript ${\rm reco}$ .

Since we have observed that fitting the field structure and $X_{\rm max}$ at the same time sometimes results in ill-converging fit, we have taken the following approach. We perform 20 reconstructions for fixed values for $X_{\rm max}$ that vary in steps of 20 g/cm² between 500 g/cm² and 900 g/cm². For each construction, we fit the three parameters of the electric-field structure $h,E$ and $\alpha$ while keeping $X_{\rm max}$ fixed. As discussed in Ref. [8], from the Stokes parameters, we know the power value and the orientation of the polarization vector, however not the sign of the electric field vector. As a result, there are two solutions that give almost the same value of $\chi^{2}$ but differ in the sign of the induced current and thus that of the electric-field direction. For this reason, we change the obtained angle $\alpha$ by 180^∘ and fit again for $h,E$ and $\alpha$ while keeping $X_{\rm max}$ fixed. We select the electric field structure and the value of $X_{\rm max}$ that give the smallest value for $\chi^{2}$ .

In Fig. 2 the difference $\Delta h=h_{\rm reco}-h_{\rm true}$ is shown as function of the true height where the electric field starts, $h_{\rm true}$ , chosen randomly between 3 km and 10 km. Lower values were not considered as the agreement for these cases is close to perfect. Per height-bin of 1 km the mean values and the standard deviations are given. For heights lower than 9 km a rather good agreement is obtained between the reconstructed height, $h_{\rm reco}$ , and the true height where the electric field starts, $h_{\rm true}$ with a standard deviation of less than 1 km. For heights between 9 km and 10 km, the spread in $\Delta h$ increases to reach a standard deviation of about 1.3 km although the mean difference is vanishingly small. The explanation is that at large altitudes there are not many particles, as the shower is still very young, and we thus loose sensitivity. Even though the spread in heights is increasing from 0.2 km at the height of 3.5 km to 1.3 km at 9.5 km it should be noted that the relative errors, $\Delta h/h_{\rm true}$ , stay more constant, ranging from 5% to 13%, in a similar range as for the single layer case.

The true value of the electric field $E_{\rm true}$ is also chosen randomly between 5 kV/m and 70 kV/m while the angle $\alpha$ is randomly selected from 0 to 360^∘. In Fig. 3 we plot $\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm true}$ as a function of the work of the net force $W=F.H$ . Here $F$ is the net force of the Lorentz force and the true electric force and $H$ is the thickness of the layer. For these one-layer cases, $H=h_{\rm true}$ . When the work of the net force of a layer is large, the amount of radio emission emitted from this layer is also large and thus we are more sensitive to this layer. As a result, the orientation of the electric field in this layer is determined more accurately (or $\Delta\alpha$ is small) for large values of $W$ . Fig. 3 shows that the angle of the field is very well reconstructed, with $\Delta\alpha=\alpha_{\rm reco}-\alpha_{\rm true}<4^{\circ}$ . For work less than 100 MeV, the deviation of $\Delta\alpha$ is about 1^∘, while for larger work of the net force, it almost vanishes.

Fig. 4 shows the difference between the reconstructed strength of the electric field and the true value, $\Delta E=E_{\rm reco}-E_{\rm true}$ . When the electric field is less than 40 kV/m, the reconstructed values from the MGMR3D calculation agrees well with the true value, with differences less than 5 kV/m. When the field strength is large, the differences become large since the current saturates at about 50 kV/m and thus we loose the sensitivity to the field strength. However, the relative errors are less than 5% when the field strenght is weaker than 50 kV/m and vary from 8% to 11% when the field strength is large.

Fig. 5 shows the difference between the reconstructed values of $X_{\rm max}$ and the true values, $\Delta X_{\rm max}=X_{\rm max}$ _reco $-X_{\rm max}$ _true. For heights lower than 7 km, there seems to be no sizable systematic differences. The mean discrepancies are less than 10 g/cm², comparable to the fair-weather results for $X_{\rm max}$ . For heights from 7 km to 10 km, the mean difference is a bit larger, from 30 g/cm² to 40 g/cm², and increases with the true height, but the standard deviation is smaller.

One other important parameter is the normalization factor $f_{r}$ for the Stokes parameters between MGMR3D and CoREAS which is given in Eq. (16). For fair-weather showers, the mean normalization factor is equal to 1. For thunderstorm showers, the mean value of the normalization factor of these 94 showers is 0.63 and the standard deviation is 0.24. This means that the radiation in MGMR3D is overestimated as compared to the CoREAS result. This could be a reflection of the fact that the reconstructed electric fields are a bit over estimated, as shown in Fig. 4.

As was shown in this section, there are no sizable systematic differences between the reconstructed parameters for the electric field structure and $X_{\rm max}$ using a MGMR3D fit and the Monte-Carlo truth as expressed using CoREAS. The systematic differences are mostly smaller than the reconstruction accuracy.

III.2.3 Multiple layers

To obtain test results for semi-realistic cases we have simulated showers passing through a two-layer electric field. Layer 1 extents from height $h_{1}$ to the height $h_{2}<h_{1}$ . Layer 2 extents from height $h_{2}$ to the ground. The height for layer 1, $h_{1}$ , is selected in the range from 5 km to 10 km. The height for layer 2, $h_{2}$ , is chosen in the range from 2 km to 4.5 km. In reality the field strength in layer 1 is usually larger than the one in layer 2. Thus, the field strength in layer 1, $E_{1}$ , and layer 2, $E_{2}$ , are chosen randomly from 2 kV/m to 60 kV/m and from 2 kV/m to 30 kV/m, respectively, each with arbitrary angles $\alpha_{1}$ and $\alpha_{2}$ with respect to the ${\bf{v}\times\left(\bf{v}\times\bf{B}\right)}$ -axis. The field strength $E_{2}$ is chosen smaller than $E_{1}$ since in reality the electric field in the cloud is usually stronger than the one between the cloud to ground. Within these constraints the parameters are chosen randomly to simulate 91 vertical showers in CoREAS. In order to find the reconstructed values using MGMR3D, following a similar procedure as was taken for the one-layer case, we perform 20 separate parameter searches with different values of $X_{\rm max}$ . We fit the six parameters of the two-layer electric field while keeping $X_{\rm max}$ fixed. We then repeat the procedure changing the values of $\alpha$ by 180^∘.

Fig. 6 shows the difference between the reconstructed and the true height, $\Delta h=h_{\rm reco}-h_{\rm true}$ , for both layers. We observe very similar differences, of the same order of magnitude as shown in Fig. 2 for the single-layer case. For small heights in layer 2 the standard deviation is of the order of 100 m, but is larger, of the order of 1000 m, for larger heights in layer 1. On average the top height is reconstructed at a too low altitude. Similar to what is seen in the single-layer case, although the spread in heights for layer 1 is much larger than that for layer 2, the relative errors vary less, from 5% to 13%.

Since we have found that the difference between reconstructed and true electric field strength correlates most strongly with the upper height of the layer, we show in Fig. 7 $\Delta E=E_{\rm reco}-E_{\rm true}$ as function of the height. We observe that the standard deviation of $\Delta E$ differs for the two layers. While the mean $\Delta E$ is constant and slightly negative for the lower layer, the mean value for $\Delta E$ is positive and decreases with height for the upper layer. The spreads in layer 1 are larger than that in layer 2 because the field strength in layer 1 is stronger than the one in layer 2 and for large field strength we loose sensitivity as shown in the one-layer case.

Fig. 8 shows that there is a correlation between $\Delta\alpha$ and the work of the net force $W$ . When the work of the net force of a layer is larger than 50 MeV, which is usually seen in layer 1, the amount of radio emission emitted from this layer is large and thus we are quite sensitive to this layer. Therefore, the orientation of the electric field in this layer is determined rather accurately or $\Delta\alpha$ is small. In contrast, when the work of the net force is smaller than 50 MeV, which is usually seen in layer 2, the amount of radiation from this layer is small and thus we loose sensitivity to the polarization or the spread in $\Delta\alpha$ is large.

Fig. 9 shows that the differences for $X_{\rm max}$ are on average positive with a mean difference of about 50 g/cm², about equal to the standard deviation. The fact that $X_{\rm max}$ is deeper in the atmosphere is probably correlated with the, on average, negative values for $\Delta h$ shown in Fig. 6.

Similar to the one-layer showers, the mean norm factor and its standard deviation of the two-layer showers are $0.62$ and $0.19$ , respectively.

To conclude, in general the reconstructed values are scattered around the true values and the strongest bias is seen in the obtained values for $X_{\rm max}$ that are, however, of little interest in the study of atmospheric electric fields.

IV Extracted electric field configuration from LOFAR data using MGMR3D

To test the proposed approach for reconstructing the atmospheric electric field from data, we have applied the reconstruction procedure to LOFAR data for a cosmic ray event that has been recorded under thunderstorm conditions [8]. We selected an event for which the radio footprint is particularly complicated, to the extent we could not reach satisfactory results in an earlier application of the procedure, event 1 from Ref. [8]. The shower axis for this event is inclined with a zenith angle of $\theta=39.2^{\circ}$ . As shown in Ref. [8], there were large differences between CoREAS and the reconstruction for this event using the older version of MGMR3D.

We have applied the present reconstruction method to recover the electric field for this event. Similar to what has been done in Ref. [8], we model the electric field with three layers. The core position is kept fixed at the position estimated from the particle-detector array LORA [25]. We perform 20 reconstructions for fixed values for $X_{\rm max}$ that vary in steps of 20 g/cm² between 500 g/cm² and 900 g/cm². For each construction, we fit nine parameters of the electric-field structure while keeping $X_{\rm max}$ fixed. At the end we select the electric field structure and the value of $X_{\rm max}$ that give the best fit result. This approach is similar to the reconstruction procedure used in Section III.2.2 and Section III.2.3. The results of the MGMR3D calculation that gives the best fit is shown in Fig. 10. The electric field reconstructed from this calculation is plugged in CoREAS to confirm the procedure. Since in CoREAS, $X_{\rm max}$ is an output which cannot be chosen before running simulations, we simulate 20 showers with the reconstructed electric-field structure and select the simulation which gives the smallest value of $\chi^{2}$ , shown in Fig. 11. The values of the extracted electric field parameters, the value of $X_{\rm max}$ , the reconstructed energy and the values of $\chi^{2}$ are given in Table. 3.

With the new version of MGMR3D, the event can be reconstructed rather well by a three-layer electric field. As shown in Fig. 10 and Fig. 11, the differences between MGMR3D and CoREAS are relatively small. It can be seen from the first panels of these two figures that Stokes I, i.e. the radiation intensity, given by both MGMR3D and CoREAS fits the LOFAR data well, except for distances near the core of the shower. At these close distances, the intensity given by CoREAS is a bit lower than the data, with the peak shifted to a slightly larger distance. The intensity at small distances is rather small because of a destructive interference between the radiation from layer 2 and layer 3. The peak in the Stokes I at a distance of 175 m from the shower axis is reproduced well in the MGMR3D calculation. Stokes $Q$ and Stokes $U$ represent the linear polarization of the radiation. The second and the third panels of Fig. 10 and Fig. 11 show the ratio $Q/I$ and $U/I$ , respectively. The polarization of the radiation is complicated, changing with the distances from the shower axis. Stoke $V$ , presenting the circular polarization, is plotted in the last panels of Fig. 10 and Fig. 11 in terms of the ratio $V/I$ . There is a large amount of circular polarization which is reflected in the change in the direction of the electric field between three layers. The linear polarization and the circular polarization are also reproduced quite well except at small distances, within 100 m from the core of the shower where the radio intensity is rather small. As shown in Table. 3 the MGMR3D and CoREAS results correspond to an energy of the cosmic ray of $1\cdot 10^{8}$ GeV and $1.07\cdot 10^{8}$ GeV. This energy differs from what is deduced from the particle detectors, $2.67\cdot 10^{7}$ GeV. The reason for this could be that the response of the LORA detectors is affected by the thunderstorm.

Since we can determine the components of the electric fields perpendicular to the direction of the shower, we can derive the purely horizontal components of the field. We decompose $\textbf{E}_{\perp}$ into two components $\textbf{E}_{\textbf{vxz}}$ and $\textbf{E}_{\textbf{vx(vxz)}}$ along vxz and vxz and shown in Table. 4. The purely horizontal components $\textbf{E}_{\textbf{vxz}}$ determined for this event are small because this shower is inclined, with the zenith angle of $\theta=39.2^{\circ}$ , since one expects the atmospheric electric field to have a strong vertical component. This component in layer 1 is larger than that in other lower layers, which is one would expect because the charge layers are not purely horizontal or the event occurred at the edge of the charged layer. The heights where the electric fields change could be the positions of different charged layers. There could be a negative-charge layer at 3.38 km and positive-charge layers at 8 km and 1.85 km altitude. This is similar to what has been determined from Lightning Mapping Array (LMA) observations of a flash that occurred near the LOFAR core in June 2019. Details about the interpretation of this flash can be found in Ref. [8].

The result of performing a reconstruction for this complicated-footprint shows that we can get a stable result from the new version of MGMR3D that captures the main structures seen in intensity and polarization. Using the reconstructed fields in the microscopic calculation gives results for the complicated radio footprint that show the same, very non-trivial, structures.

Layer	h [km]	E [kV/m]	$\alpha$ [^∘]
1	8.012	34.900	157.3
2	3.381	99.377	-64.6
3	1.851	71.756	102.2
$X_{\rm max}$ _3D[g/cm²]			860
$X_{\rm max}$ _C[g/cm²]			848
Energy_3D[GeV]			$1.0\times 10^{8}$
Energy_C[GeV]			$1.07\times 10^{8}$
$\chi^{2}_{\rm 3D}$			$2.39$
$\chi^{2}_{\rm C}$			$3.76$

Table 3: Nine extracted electric field parameters,

X_{\rm max}

and the reconstructed energy of event 1

Layer	h [km]	$\textbf{E}_{\textbf{vxz}}$ [kV/m]	$\textbf{E}_{\textbf{vx(vxz)}}$ [kV/m]
1	8.012	24.721	24.635
2	3.381	-5.547	-99.222
3	1.851	-12.460	70.665

Table 4: The components

\textbf{E}_{\textbf{vxz}}

and

\textbf{E}_{\textbf{vx(vxz)}}

of the electric fields determined from event 1

V Summary and Conclusions

In this work, we have made an comprehensive comparison between the radio emission during thunderstorm conditions from MGMR3D and from CoREAS. As a result, we improve the parameterization for MGMR3D in the presence of atmospheric electric fields.

With the parameters determined, we have generated a large number of radio footprints for showers passing through a layered electric field using the microscopic code CoREAS. From these footprints we have reconstructed the electric field structure using MGMR3D and compared these with the true structure that was used for generating the footprint. This shows that the layer heights can be reconstructed with an accuracy of better that $\pm 1$ km, the field strengths with $\pm 10$ kV/m, and orientation angles within a few degrees. Exceptions are the orientation of weak fields for cases in which there are also strong fields involved as well as cases where the top height of the electric field layer lies at a height where the shower is still young with relatively few charged particles. In the region where we are sensitive to the parameters of the electric fields, the differences between MGMR3D and CoREAS are rather small for all vertical showers passing through a one-layer electric field or a two-layer one. Large discrepancies are observed for the cases where we loose sensitivity to the parameters of the electric fields seen in both MGMR3D and CoREAS.

We have applied the approach to extract electric field structure to an event having a complicated radio footprint as measured by LOFAR and we are able to reconstruct the main features of this event. Therefore, it can be concluded that MGMR3D can be used to reconstruct the structure of the electric fields by using the radio emission emitted from air showers passing through thunderclouds. This method of determining atmospheric electric fields will help to study about the process of lightning initiation and propagation.

Acknowledgement

This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.01-2019.378. K. D. de Vries was supported by the European Research Council under the EU-ropean Unions Horizon 2020 research and innovation programme (KdV, grant agreement No 805486).

Appendix A Parameter values

Parameter	$J_{Q}$	$X_{0}$	$\lambda$
Value	0.21	100	100

Table 5: Fixed values of parameters

Appendix B Programming details

The program can be downloaded from MGMR3D-v3 as a zip file. Make sure you run version 3.

References

Schellart et al. [2015] P. Schellart et al., Probing atmospheric electric fields in thunderstorms through radio emission from cosmic-ray-induced air showers, Phys. Rev. Lett. 114, 165001 (2015).
Scholten et al. [2008] O. Scholten, K. Werner, and F. Rusydi, A macroscopic description of coherent geo-magnetic radiation from cosmic-ray air showers, Astroparticle Physics 29, 94 (2008).
Werner et al. [2012] K. Werner, K. D. de Vries, and O. Scholten, A realistic treatment of geomagnetic cherenkov radiation from cosmic ray air showers, Astroparticle Physics 37, 5 (2012).
Scholten and Werner [2009] O. Scholten and K. Werner, Macroscopic model of geomagnetic-radiation from air showers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 604, S24 (2009), {ARENA} 2008.
Askaryan [1962] G. Askaryan, Excess Negative Charge of the Electron-Photon Shower and Coherent Radiation Originating from It. Radio Recording of Showers under the Ground and on the Moon, J. Phys. Soc. Japan Vol: 17, Suppl. A-III (1962).
de Vries et al. [2010] K. D. de Vries et al., The lateral distribution function of coherent radio emission from extensive air showers: Determining the chemical composition of cosmic rays, Astroparticle Physics 34, 267 (2010).
Trinh et al. [2016] T. N. G. Trinh et al., Influence of atmospheric electric fields on the radio emission from extensive air showers, Phys. Rev. D 93, 023003 (2016).
Trinh et al. [2020] T. N. G. Trinh, O. Scholten, S. Buitink, U. Ebert, B. M. Hare, P. R. Krehbiel, H. Leijnse, A. Bonardi, A. Corstanje, H. Falcke, T. Huege, J. R. Hörandel, G. K. Krampah, P. Mitra, K. Mulrey, A. Nelles, H. Pandya, J. P. Rachen, L. Rossetto, C. Rutjes, S. ter Veen, and T. Winchen, Determining electric fields in thunderclouds with the radiotelescope lofar, Journal of Geophysical Research: Atmospheres 125, e2019JD031433 (2020), e2019JD031433 10.1029/2019JD031433.
Alvarez-Muñiz et al. [2012] J. Alvarez-Muñiz, W. R. C. Jr., and E. Zas, Monte carlo simulations of radio pulses in atmospheric showers using {ZHAireS}, Astroparticle Physics 35, 325 (2012).
Huege et al. [2013] T. Huege, M. Ludwig, and C. W. James, Simulating radio emission from air showers with coreas, AIP Conference Proceedings 1535, 128 (2013).
Zilles et al. [2020] A. Zilles, O. Martineau-Huynh, K. Kotera, M. Tueros, K. de Vries, W. Carvalho, V. Niess, N. Renault-Tinacci, and V. Decoene, Radio Morphing: towards a fast computation of the radio signal from air showers, Astropart.Phys. 114, 10 (2020).
Marin and Revenu [2012] V. Marin and B. Revenu, Simulation of radio emission from cosmic ray air shower with {SELFAS2}, Astroparticle Physics 35, 733 (2012).
Scholten et al. [2018] O. Scholten, T. N. G. Trinh, K. D. de Vries, and B. M. Hare, Analytic calculation of radio emission from parametrized extensive air showers: A tool to extract shower parameters, Phys. Rev. D 97, 023005 (2018).
Schellart et al. [2014] P. Schellart et al., Polarized radio emission from extensive air showers measured with lofar, Journal of Cosmology and Astroparticle Physics 2014 (10), 014.
Bossard et al. [2001] G. Bossard, H. J. Drescher, N. N. Kalmykov, S. Ostapchenko, A. I. Pavlov, T. Pierog, E. A. Vishnevskaya, and K. Werner, Cosmic ray air shower characteristics in the framework of the parton-based gribov-regge model nexus, Phys. Rev. D 63, 054030 (2001).
Mitra [2021] P. Mitra, High Precision Reconstruction of Air Shower Properties With Dense Radio Arrays, Ph.D. thesis, Vrije Universiteit Brussel (2021).
Heck et al. [1998] D. Heck et al., CORSIKA: a Monte Carlo code to simulate extensive air showers (1998).
S. Buitink et al. [2007] S. Buitink et al., Amplified radio emission from cosmic ray air showers in thunderstorms, A&A 467, 385 (2007).
Buitink et al. [2014] S. Buitink et al., Method for high precision reconstruction of air shower ${X}_{\mathrm{max}}$ using two-dimensional radio intensity profiles, Phys. Rev. D 90, 082003 (2014).
van Haarlem et al. [2013] M. P. van Haarlem et al., Lofar: The low-frequency array, A&A 556, A2 (2013).
[21] Terminal velocity, https://www.grc.nasa.gov/WWW/K-12/airplane/termv.html.
Gaisser and Hillas [1977] T. K. Gaisser and A. M. Hillas, Reliability of the Method of Constant Intensity Cuts for Reconstructing the Average Development of Vertical Showers, in International Cosmic Ray Conference, International Cosmic Ray Conference, Vol. 8 (1977) p. 353.
Salvador et al. [2021] A. Salvador, N. Pineda, J. Montanyà, J. A. López, and G. Solà, Thunderstorm charge structures favouring cloud-to-ground lightning, Atmospheric Research 257, 105577 (2021).
Dwyer and Uman [2014] J. R. Dwyer and M. A. Uman, The physics of lightning, Physics Reports 534, 147 (2014), the Physics of Lightning.
Thoudam et al. [2014] S. Thoudam et al., Lora: A scintillator array for {LOFAR} to measure extensive air showers, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 767, 339 (2014).