Analytic Calculation of Radio Emission from Extensive Air Showers subjected to Atmospheric Electric Fields

We have developed a code that semi-analytically calculates the radio footprint (intensity and polarization) of an extensive air shower subject to atmospheric electric fields. This can be used to reconstruct the height dependence of atmospheric electric field from the measured radio footprint. The various parameterizations of the spatial extent of the induced currents are based on the results of Monte-Carlo shower simulations. The calculated radio footprints agree well with microscopic CoREAS simulations.


Introduction and summary
Atmospheric electric fields during thunderstorm conditions can drastically modify the radio footprint of extensive air showers (EAS). This fact has been used [1,2] to extract from the measured radio footprint the altitude dependence of the electric fields. In contrast to conventional weather-balloon measurements this offers a non-intrusive way to determine these fields that are essential for understanding lightning initiation and development [3].
The extraction of the atmospheric fields relies on a model calculation of radio emission from an EAS. Accurate codes exists that perform reliable microscopic calculations [4]. One such code that allows for setting arbitrary atmospheric electric fields is CoREAS [5] that uses results of a CORSIKA shower evolution. This uses a Monte-Carlo simulation and thus the results of two calculations with very similar but unequal electric fields may differ considerably due to random shower-to-shower fluctuations. As a result it is close to impossible to use CoREAS in a systematic optimization procedure. To resolve this problem we have developed a deterministic semi-analytic calculation that predicts the radio-footprint given a certain altitude dependence of the induced currents in the EAS. An additional advantage is that the analytic code is much less CPU intensive than the microscopic calculation (a few seconds compared to many hours). As a result it now becomes practical to solve the inverse problem from radio emission, i.e. to determine from a measured radio footprint (intensity and polarization distributions on the ground) of a shower the configuration of currents in the atmosphere that produced this.
The code is based on the macroscopic formulation of radio emission as developed in refs [6][7][8][9][10] and incorporated in the EVA-code [11]. An important difference between the EVA-code and the e-mail: Scholten@KVI.nl present code (called EVA-light) is that here we use a direct parametrization of the current distribution in EAS. Certain parameters, such as the distribution of the current density as function of distance to the shower axis, are kept fixed, while other parameters related to the altitude dependence of the currents are optimized to fit the data. In addition the new code is optimized for a fast calculation of the complete radio footprint.
Work is still ongoing to refine the parametrizations of the charge-excess and the induced-current clouds where the realistic distributions as determined in Monte-Carlo shower evolution calculations are guiding. Preliminary results show good agreement with those of microscopic CoREAS calculations.

Modelling Radio emission from EAS
The currents and charges in the EAS are modeled as a charge cloud with a parameterized density profile moving towards Earth with the light velocity. Radio emission is subsequently calculated using the retarded vector potential where we follow the a macroscopic approach to construct the retarded potentials. Furthermore we have put some effort in optimizing the code for the calculation of a complete radio footprint.
In this paper we work with the four-current j µ (t s , x s , y s , h) where z s = −t s (taking c = 1) is the distance from the shower front to the ground, called height and always measured along the shower axis, x s , and y s are the transverse directions, and h is the distance from the shower front to the point in the charge cloud. We use the notation where µ = 0 denotes the time (charge) components and µ = x, y, z the space (current) components of a four vector. The shower front is moving with the light velocity. A particular point in the charge cloud is thus at an height of ζ = z s + h.
Following the usual notation, the vector potential for an observer at a point (x o , y o , z o = 0) on the shower plane, defined as the plane perpendicular to the shower axis going through the point of impact of the shower on the ground, is taken as where L is the optical path-length, L = c(t o − t r ). Only for a homogeneous medium with a constant index of refraction n it is given by for a moving point charge with velocity cβ. The vector potential can now be written more compactly as The index of refraction depends on the height in the atmosphere through it density dependence as given by the Gladstone-Dale relation, n GD = 1 + n ρ ρ(z) where ρ is the air density. For the evaluation of the retarded distance given in Eq. (2), the average value of the index of refraction over the photon trajectory will be used assuming the straight-line approximation for the photon path [11], which does not depend on the distance of the observer to the shower axis.
To evaluate the integration over z in Eq. (3), we follow the same approach as used in Ref. [11] and replace the integral in the z-direction by an integral over λ = h c − h where D = 0 at the critical height h c . This substitution allows for an easier calculation of derivatives of the vector potential that are necessary to calculate the electric field since the derivatives vanish on the limits of the λ integration. As shown in Ref. [11]

Electric field, charge excess
The charge excess in the shower front is given by Q(z) propagating with the light velocity in the −ẑ-direction thus contributing to the zero and the z component of the vector potential. Substituting the expression for the vector potential, integrating over the spatial extent of the charge cloud and converting to the integration over λ, the radial component of the electric field can be written as where r 2 s = (x 2 s + y 2 s ) is an integration point in the shower front and the position of the observer is denoted with the subscript o. The integral is separated into an integration over λ, a line-sub-shower, along a line parallel to the axis of the full shower, followed by one over the radial extent. In the numerical calculation the results line-sub-showers, I CX λ (t), are stored for a grid of r s and d values, where the latter is the distance of the observer to line-sub-shower.
In order to obey charge conservation a residue charge line-density Q (z) = (dQ/dz) is left behind at rest in the atmosphere which can be ignored as shown in Ref. [7].

Electric field, Transverse current
Along a similar line of reasoning as for the charge-excess radiation, the transverse current contribution is written as where J x = dI x /dt r and f (h) = d f /dh. As before the results for line-sub-showers, I TC λ (t) are calculated once on a grid of r s and d values to be used subsequently in the calculation of the complete footprint of the electric field.
Due to the continuity equation we have dJ x /dx = dQ/dt, thus for r 2 and a similar term for J y . This corresponds to a di-pole charge distribution moving with the shower front. In reality these particles thermalize and will trail well behind the front and thus not contribute to coherent emission. To account for this we include a decrease proportional to e −X/a D X 0 , where X depends on t through its dependence on the height in the atmosphere. In Eq. (8) X 0 is the electron mean free path and a D a parameter that is adjusted to get best agreement with microscopic shower calculations. Following a similar notation as before the contribution to the electric field is ; where a similar expression is obtained for E D y .

Parameterizations
The density distribution of the four current is parameterized as (10) where r = x 2 + y 2 , u = h/α λ, M 0 the Moliere radius, h the distance behind the shower front. The normalization constants N w and N f in functions w and f are determined by w(r) dr = 1 and ∞ 0 f (h, r) dh = 1 ∀ r. In this way J µ (t) is the charge and current at time t integrated over the whole shower front. Note that w(r) corresponds to r times the NKG function for s = 2. The two parameters in f , the pancake thickness λ(r) and α(X), depend on distance to the shower axis, r, and penetration depth, X, as λ(r) = max[λ 0 , λ 1 (r/r 1 )] ; α(X) = 1 + 0.5 × X/X max .
The radial dependence of the pancake thickness is given by λ(r) such that for r = 0 we have λ(0) = λ 0 increasing linearly with distance on an exponential grid where at a distance of r 1 =100 m from the core we have λ(100) = λ 1 . With increasing α the maximum in the particle density shifts to larger distances behind the shower front. Optimization of the parameters is still ongoing. We note that choosing M 0 = 30 m, λ 0 = 0.15 m, λ 1 = 3 m, reproduces the main features of the charge cloud as deduced from microscopic calculation. This also gives rise to a radio footprint that agrees well with that obtained from CoREAS simulations. We are still investigating the extent to which the distributions need to be chosen differently for charge excess and transverse current in view of our findings in Ref. [2]. We note that changing the pancake thickness with distance to the shower axis is essential to be able to reproduce the (not frequency filtered) pulse shape as obtained from the CoREAS calculations. The penetration-depth dependence of α does not have a strong influence on the calculations.