Analytic calculation of radio emission from parameterized extensive air showers, a tool to extract shower parameters

The radio intensity and polarization footprint of a cosmic-ray induced extensive air shower is determined by the time-dependent structure of the current distribution residing in the plasma cloud at the shower front. In turn, the time dependence of the integrated charge-current distribution in the plasma cloud, the longitudinal shower structure, is determined by interesting physics which one would like to extract such as the location and multiplicity of the primary cosmic-ray collision or the values of electric fields in the atmosphere during thunderstorms. To extract the structure of a shower from its footprint requires solving a complicated inverse problem. For this purpose we have developed a code that semi-analytically calculates the radio footprint of an extensive air shower given an arbitrary longitudinal structure. This code can be used in a optimization procedure to extract the optimal longitudinal shower structure given a radio footprint. On the basis of air-shower universality we propose a simple parametrization of the structure of the plasma cloud. This parametrization is based on the results of Monte-Carlo shower simulations. Deriving the parametrization also teaches which aspects of the plasma cloud are important for understanding the features seen in the radio-emission footprint. The calculated radio footprints are compared with microscopic CoREAS simulations.


I. INTRODUCTION
When a high-energy cosmic particle impinges on the atmosphere of Earth, it creates an extensive air shower (EAS). The electrons and positrons in the plasma cloud at the shower front are deflected in opposite directions due to the Lorentz force caused by the geomagnetic field. This creates an electric current. Since the number of particles in the EAS changes with penetration depth, the electric current in the plasma cloud changes as a function of height in the atmosphere. This varying current emits radio waves [1][2][3][4] where the intensity pattern on the ground, the intensity footprint, depends on the variation of the current with height. The penetration depth where the particle number reaches its maximum, X max , strongly depends on the specifics of the first collisions, in particular their multiplicity and thus the mass of the initiating cosmic ray. Different values of X max result in differences in the longitudinal variation of the currents which is reflected in the intensity footprint. Thus X max can be reconstructed on the basis of the footprint which allows for a determination of the mass composition of cosmic rays [5] for fair-weather conditions.
In Refs. [6,7], a new method is introduced to determine the electric fields in the atmosphere during thunderstorm conditions by using the radio footprint from an EAS. The basic principle is the same as used in determining X max for air showers recorded under fair-weather conditions (fairweather events). During thunderstorm conditions an additional strong variation of the current with height is induced by the thunderstorm electric field [6,7], which also varies with height in direction and magnitude. During such conditions the effect of the atmospheric electric field on the current can be dominant. While it is sufficient to consider only the intensity footprint for fair-weather showers, the footprints for all Stokes parameters [8,9] are necessary for a complete mapping of the fields in the atmosphere [10,11].
Only a single parameter, X max , needs to be extracted from the radio footprint for a fair-weather event. However, for a shower recorded under thunderstorm conditions (thunderstorm event) there are many more, order of 10, where the precise number of parameters depends on the level of sophistication of the modeling of the electric fields in the atmosphere. Therefore a simple grid search algorithm suffices to extract the value of X max for a fair-weather event, while such a grid search is totally impractical for a thunderstorm event. To make such a parameter search more efficient, one needs to be able to deterministically calculate the radio footprint given the structure of the shower, such that an infinitesimal change in the shower parameters results in an infinitesimal change in the radio footprint. In addition, it is convenient if a single calculation takes little computer resources, as such a calculation has to be iterated many times to find the optimum. These two conditions are not met by the presently available microscopic codes, CoREAS [12] and ZHAireS [13]. Since both of these codes are based on a Monte Carlo calculation of the EAS, changing a single shower parameter will, in general, affect the complete shower evolution in a nondeterministic way. Furthermore, the codes are rather computer intensive as they work on a microscopic level, tracing the individual electrons.
One approach [14] to this inverse problem is to use a CoREAS calculation as a template from which the emission amplitudes from the different shower slices is stored. The emission from other showers is calculated from the template by simply adjusting the height-dependence of the weighting factor of the shower slices extracted from the template shower. Although very promising, several details still need further attention before this procedure can be applied.
As an alternative approach to solving the posed inverse problem we have developed a code that semianalytically calculates the radio footprint of an extensive air shower for an arbitrary longitudinal structure of the electric current density in the shower front. The analytic calculation uses a negligible amount of computer time, it is about 4 orders of magnitude faster [15] than a microscopic calculation (at E ¼ 10 16 eV), and, most importantly, does not suffer from random showerto-shower fluctuations. Therefore it can be used in a chisquare minimization procedure to obtain the longitudinal structure that best reproduces the measured footprint.
For the present analytic code constructions similar to those that have been developed for the EVA code [3] are used to obtain the radiation fields from the Liénard-Wiechert potentials. Like EVA, the present code accounts for the proper retardation effects. In EVA the parametrization of the plasma cloud is obtained from a Monte Carlo simulation of a shower to be able to account for shower-toshower fluctuations to the full extent. In the present code, however, we want to eliminate all shower-to-shower fluctuations, and the shower evolution, including the structure of the plasma cloud, is parametrized. In this respect the model is similar to the MGMR model [1], with the important difference that here the radial extent of the plasma cloud is taken into account. For this reason we named it MGMR3D. Charge-excess radiation and a realistic index of refraction of air are also taken into account.
It is known from shower universality that the largest shower-to-shower fluctuations occur in the longitudinal shower evolution whereas the structure of the plasma cloud, like lateral extension and thickness, shows hardly any shower-to-shower fluctuations. In order to keep the number of parameters manageable, we have, therefore, adopted the approach where we use a generic parametrization of the structure of the plasma cloud where the dependence of the plasma cloud integrated currents on the height in the atmosphere is left free since this is what we want to extract from the radio footprint. The parametrization of the plasma cloud is inspired by the results of Monte Carlo shower simulations and discussed in Sec. II. To validate the parametrizations, we have verified in Sec. III that the results of the MGMR3D calculations for the radio footprint agree sufficiently well with those from microscopic CoREAS simulations. Obtaining such a parametrization is an important part of the present work. As an interesting spin-off, the code also allows one to investigate which are the essential parameters of the plasma cloud for certain features of the radio pulse footprint. In this way-as an example-we noted that the temporal structure of the radio pulse, a strong peak followed by a very shallow undershoot, is strongly determined by the radial dependence of the pancake thickness, see Sec. III.
Since the present code, MGMR3D, is supposed to facilitate an iterative approach to reproduce a measured radio footprint, much attention was given to its numerical stability and its calculational speed. The code, MGMR3D, is available from the authors upon request.

II. MODELING RADIO EMISSION FROM EAS
The currents and charges in the EAS are modeled as a plasma cloud with a parametrized density profile moving towards Earth at the speed of light, c. These currents are used to construct the retarded Liénard-Wiechert potential from which the radiation fields are calculated. Here we closely follow the approach used in the EVA code.
To parametrize the charge cloud we introduce the shower-fixed coordinates ðt s ; x s ; y s ; hÞ where t s is the time when the shower front is at a distance z s ¼ −t s c from the ground (measured along the shower axis), and ðx s ; y s Þ are transverse coordinates of a point in the shower plasma cloud at a distance h behind the shower front. The structure of the charge and current distributions are expressed through a four-current j μ ðt s ; x s ; y s ; hÞ. We use the notation where μ ¼ 0 denotes the time (charge) components and μ ¼ x, y, z the space (current) components of a four vector. A particular point in the charge cloud is at a height of ζ ¼ z s þ h as measured along the shower axis.
Following the usual notation where t r denotes retarded time, the vector potential for an observer at a point ðt o ; x o ; y o ; z o Þ in the shower plane, defined as the plane perpendicular to the shower axis going through the point of impact of the shower on the ground (z o ¼ 0), is taken as Here L is the optical path length, L ¼ cðt o − t r Þ, which for a homogeneous medium with a constant index of refraction n is given by L ¼ nR ¼ nj ⃗ x o − ⃗x 0 j. Following the approach used in EVA [3], we introduce D ¼ Ljdt o =dt r j as the retarded distance. The vector potential can now be written more compactly as In the Appendix, Secs. A 1 and A 2, more details are given on the numerical calculation of the radiation fields from the vector potential Eq. (2). For simplicity the dependence of the index of refraction on the height in the atmosphere is taken as given by the Gladstone-Dale [16] relation, where n ρ is chosen such that the refractivity equals dn ¼ n − 1 ¼ 0.0003 at ground level. This can be replaced by a more realistic dependence such as used in Ref. [17] that takes the dependence of the refractivity on temperature and air humidity into account. For the evaluation of the retarded distance the average value of the index of refraction will be used assuming the straight-line approximation for the photon path [3], which does not depend on the distance of the observer to the shower axis. Thus, to a good approximation (1:straight line trajectory; 2:index of refraction does not change with observer position) the retarded distance can be expressed as for a moving point charge with velocity βc and an observer at time t o at a distance d ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The four current is parametrized as j μ ðt s ; x s ; y s ; hÞ ¼ wðr s Þ r s fðh; r s ÞJ μ ðt s Þ; where r s ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi x 2 s þ y 2 s p , and the functions w and f are normalized according to R wðrÞdr ¼ 1 and R ∞ 0 fðh; rÞdh ¼ 1 ∀ r. In this way J μ ðt s Þ is the charge and current at time t s integrated over the complete plasma cloud of the EAS at a fixed time.

A. Parametrizations
In the parametrizations of the charge cloud we first concentrate on the structure of the charge/current cloud under fair-weather circumstances. The effect of atmospheric electric fields is considered in a following step. It should be noted that these fields will primarily change the magnitude and orientation of the transverse electric currents, J μ ðt r Þ; however, as shown in Ref. [7], these fields will also affect the structure of the current cloud, in particular fðhÞ, see Eq. (6).
Most of the plasma-cloud parameters are obtained through a comparison with the results of a Monte Carlo simulation of an extensive air shower. Based on shower universality we use a single shower in order not to have complications of averaging the results of two showers with different values for X max . As noted in the following discussions, for some cases we have preferred a CONEX-MC simulation, the same that lies at the basis of the EVA calculations [3], for ease of extracting more detailed information of the shower structure while for others, where the atmospheric electric field is important, we used CORSIKA.

Fair-weather conditions
The spatial extent of the charge cloud is modeled as given in Eq. (6) where the functions w and f are normalized to unity. Under fair-weather conditions the radial dependence of the transverse current is parametrized as with ξ ¼ r s =M 0 introducing the Moliere radius M 0 as a scaling factor and where N w is chosen such that R wðrÞdr ¼ 1. Note that wðr s Þ corresponds to r s times the Nishimura-Kamata-Greisen function for s ¼ 2 [18]. In Fig. 1 we compare the results of CONEX-MC [19] (the Monte Carlo version of CONEX) with the parametrization Eq. (7) using the parameters given in Appendix B.
The current density at a distance h behind the shower front is parametrized as where η ¼ h=λ. The norm, N f is chosen such that R ∞ 0 fðh; r s Þdh ¼ 1 for all values of r s . The pancake-thickness scaling parameter, is factorized in a dependence on distance to the shower axis, r s , and a scale parameter αðjFjÞ that depends on the net force acting on the electrons, jFj, and is specified in more detail in Eq. (19). The radial dependence of the pancake thickness is such that near the shower core we have Λð0Þ ¼ Λ 0 while increasing linearly at larger distance, where at a distance of In principle the pancake-thickness scale parameter α may also depend on penetration depth. We observe however that for fair-weather showers the radio footprint as well as pulse shape (discussed in the following chapter) show very little sensitivity to such a height dependence and are determined almost solely by the pancake thickness at the shower maximum. We thus ignore this possible dependence.
In Fig. 2 the parametrization of the pancake structure, using parameter values from Appendix B, is compared to the results of a CONEX-MC [19] simulation for a fair-weather shower at the shower maximum.
The parametrization of the longitudinal shower profile for the charge excess and the transverse current in the plasma cloud for a shower under fair-weather conditions is based on the Gaisser-Hillas formula [20] for the dependence of the number of charged particles on X z , the penetration depth, where γ is a parameter controlling the width of the distribution and X 0 the reference point. The transverse current, see Eq. (6), is obtained by multiplying the number of charged particles with the drift velocity, where the induced transverse drift velocity is denoted as ⃗u ⊥ .
The drift velocity will increase with increasing forces acting on the charges; however, for large forces, as we will encounter during thunderstorm conditions, one should be careful to take into account that the velocity of the particles does not exceed the speed of light. Following the arguments given in Ref. [7] we take, where the parameter υ 0 , as discussed in a later section, is taken such that for fair-weather conditions one is still in the regime where the drift velocity scales linearly with the Lorentz force. When not too large, the drift velocity is proportional to the force acting on the plasma charges, where F β represents a friction constant, a t ¼ 2, and a normalization constant of X v ¼ 500 g=cm 2 is used. The last factor in Eq. (14) takes into account the fact that the drift velocity depends on the penetration depth in the atmosphere. At low altitudes the drift velocity decreases due to increased density and at high altitudes, early in the  shower development, the energy of the particles in the plasma is enormous, and for this reason their sidewards drift is small. For fair-weather conditions only the Lorentz force is acting, ⃗ F ⊥ ¼ e ⃗v s × ⃗ B, and is constant along the shower. Here ⃗v s is along the shower with magnitude c, and ⃗ B is the magnetic field. In Fig. 3 the parametrization for the drift velocity is compared with the results of a CONEX-MC calculation for a vertical shower where a geomagnetic field of 56 μT is used at an angle of 63°with the horizontal. Using this figure the values for F β and a t in Eq. (14) are deduced where their values are given in Appendix B. The arguments for introducing υ 0 are discussed in a following section. The longitudinal current due to charge excess is defined as, where, models the dependence of the charge excess fraction on penetration depth with a c ¼ 0.5. This parametrization is compared with the results of a simulation in Fig. 3, while J 0 Q is a normalization constant. For an inclined shower at an angle θ s , ignoring the curvature of Earth, we have H ¼ z cos θ s , and the height dependence of the atmospheric penetration depth is taken as where H denotes the height above ground, and where the constants a, b, c depend on height as given in Table I for the U.S. standard atmosphere [21], which are the same as used in CORSIKA [22].

Thunderstorm conditions
In the presence of thunderclouds the air shower will generally evolve through areas in the atmosphere where there are large electric fields. These will significantly alter the currents in the shower front [6,7]. It is precisely these currents we want to determine from the radio footprint. In leading order the electric fields will change the magnitude and direction of the induced drift velocities, in Eq. (14). We will assume that the strength of the component of the field parallel to the shower is below the limit where secondary electron avalanches are formed. For these strengths the number of particles in the shower is not significantly affected [7] and we thus can ignore it. The component of the field perpendicular to the shower directly influences the current. An important secondary effect of the electric field is to increase the pancake thickness [7]. Since the particles in the shower front are constantly being regenerated from the more energetic ones that drive the shower, the pancake structure is affected only at those heights where the field acts. The structure of the plasma cloud is hardly determined by its structure at larger heights, and it thus can be said that there are no memory effects. This has been checked with Monte Carlo simulations.
Monte Carlo simulations show that the distance from the shower front that contains 50% of the number of charged particles near the shower maximum depends quadratically on  the transverse force, see Fig. 4. This effect we have parametrized through a dependence of the pancake-hickness scaling factor α, see Eq. (9), on the strength of the transverse force, jFj, acting on the particles in the plasma cloud, The parameter a E is adjusted to reproduce the median trailing distance behind the shower front as obtained from CORSIKA simulations, see Fig. 4. Using the parameter values given in Appendix B, a good agreement is obtained.

III. COMPARISON WITH CoREAS SIMULATIONS
Having fixed the parameters of the plasma cloud as given in Appendix B we need to verify that the used parametrization is sufficiently detailed such that the radio footprint from MGMR3D agrees with the CoREAS results to a reasonable accuracy.
As a first step we compare in Fig. 5 the unfiltered pulses separately for the transverse current and the charge excess contributions at various distances to the core with the results of a CoREAS simulation. The calculations are performed assuming fair-weather circumstances for a vertical shower with X max ¼ 540 g=cm 2 . At this point it is worthwhile to note that-in principle-there can be an additional contribution to the emitted radiation due to an induced dipole distribution that is comoving with the shower front [1]. The geometry of its radiation pattern is very similar to that of geomagnetic radiation; however, the pulse is very elongated in time. We have calculated such a contribution and seen that best agreement with the microscopic calculation is obtained by setting this contribution to zero. This implies that the net displacement of electrons and positrons in the shower front is vanishingly small. This is in line with the conclusions reached in Ref. [23]. Another interesting point is that the typical structure of the pulse with a large positive peak followed by a long negative tail is intimately linked to the radial dependence of the pancake thickness. When taking a less pronounced increase with distance [i.e. decreasing the value of Λ 1 in Eq. (10)] the negative tail gets shorter and more pronounced. For the case in which λ is independent of distance to the core a clear bipolar pulse is obtained.

A. Fair-weather footprint
We will investigate the radio footprint of an air shower using Stokes parameters since these capture the complete polarization structure of the radio pulse. Because the objective of the present work is to develop a scheme that can be used to ease the interpretation of data, we construct the Stokes parameters with the Low-Frequency Array (LOFAR) [24] cosmic-ray experiment [8] in mind. This  implies filtering the signal to the 30-80 MHz band. In terms of the sampled pulse in the polarization direction p, where the complex voltage of the i th sample is denoted as E i;p ¼ E i;p þ iÊ i;p , the Stokes parameters can be expressed as i;p is the Hilbert transform [8] of the real measured voltage E i;p . The polarization directions p are taken alongê v×B and e v×ðv×BÞ which are by construction perpendicular to the propagation direction of the photon (in very good approximation). We sum over the whole trace. The linear-polarization angle with the v × B-axis, ψ, can be calculated directly from the Stokes parameters as ψ ¼ 1 2 tan −1 ðU=QÞ. The relative amount of circular polarization is given by V=I.
A comparison of the Stokes parameters [8,10] with the results of a microscopic CoREAS simulation is presented in Fig. 6 for the case of a relatively small value, X max ¼ 540g=cm 2 and in Fig. 7 for a larger value, X max ¼ 690 g=cm 2 . A simple geometry is used with a vertical shower and a horizontal magnetic field with B ¼ 40 μT There is a good agreement in the polarization directions between MGMR3D and CoREAS.
In Fig. 8 we present the results for a shower with a zenith angle of 30°and X max ¼ 693 g=cm 2 . The magnetic field is at an angle of α vB ¼ 60 degrees with the shower axis. In such a configuration care should be taken with the subtle effect that there is an additional component in the Lorentz force due to the drift velocity of the particles. This is taken into account by adding to the currents an additional component proportional to the component of the magnetic field parallel to the shower axis, ⃗ B ∥ , As can be seen from Fig. 8 the intensity footprint calculated with MGMR3D agrees rather well with the one calculated with CoREAS although there are some systematic differences. The reasons for these differences is not understood. Another surprising, and not understood (small) discrepancy between CoREAS and MGMR3D is seen for V=I. Even with the correction from Eq. (21) the circular polarization near the core is vanishingly small in MGMR3D while CoREAS shows a clear offset.

B. Footprint at higher frequencies
The comparisons of MGMR3D with CoREAS in the previous section have been done for the LOFAR frequency bandwidth of 30-80 MHz. To test the suitability of the macroscopic approach in other bandwidths we have calculated the radio footprint for the shower used for Fig. 6 for two other bands.
In the 100-200 MHz band, one sees from Fig. 9 that the Cherenkov ring starts to emerge with a peak intensity in a broad ring at a distance of 50-100 m. The intensity depends on azimuth angle due to the interference with charge-excess radiation. The results of MGMR3D agree quite well with those of the CoREAS simulation. At distances exceeding 200 m the agreement worsens due to the fact that the CoREAS simulation creates a noisy signal. Since this noise is independent of distance it starts to dominate over the signal at about 200 m for this calculation. As a result the relation for the Stokes parameters, Q 2 þ U 2 þ V 2 ¼ I 2 is no longer obeyed. This shows as a drop in Q=I while U=I and V=I remain small. The MGMR3D does not suffer from stochastic noise.
The Cherenkov ring is fully developed in the 200-500 MHz band as can be seen from Fig. 10. The agreement between MGMR3D and the CoREAS simulation is still very convincing, even though the differences in calculated intensities show a stronger systematic trend. At the higher frequency the problem with numerical noise is enhanced, which is why the CoREAS simulation is no longer reliable beyond 150 m distance from the axis as well as at distances less than 50 m. At a distance near 100 m the circular polarization V=I is negligible in the MGMR3D calculation while CoREAS shows considerably larger values. We have not explored the source of this difference that seems to point to an underestimate of the difference in emission heights between charge excess and transverse current radiation in MGMR3D.

C. Footprint under thunderstorm conditions
With increasing force on the electrons, the power of the emitted pulse increases until a maximum is reached for a field of the order of 50 kV=m [7]. The reason for this saturation is twofold: (a) the transverse drift velocity is limited because the velocity of the particles cannot exceed c, as expressed in Eq. (13); and (b) the pancake thickness increases with increasing field strength as shown in Fig. 4. It should be noted that these two effects are related since as the transverse velocity increases, the longitudinal component must decrease and thus the particles lag further behind the shower front. In Fig. 11 the maximal peak power  calculated using the semianalytic approach is compared with the results of a microscopic calculation. The same two-layer field configuration is used as in Ref. [7], a top layer between 8 and 3 km with a net force strength F and a lower layer from 3 km till the ground with a strength 0.3 × F in the opposite direction where the strength of F is varied. For the shower X max ¼ 500 g=cm 2 is taken. With the parameter υ 0 set to the value given in Appendix B an excellent agreement is obtained.
In Figs. 12 and 13 the footprint is compared to the results of a CoREAS simulation for two different realistic atmospheric electric field configurations as can be expected during thunderstorms as given in Table II. In Fig. 12 the electric fields in the two layers have opposite directions and the emission of the top and bottom layer interfere destructively near the core which becomes less efficient with increasing distance due to the finite refractivity of air. As a result a ringlike structure is seen in the intensity. Beyond the ring the intensity falls off a bit faster in the MGMR3D calculation as compared to CoREAS. We have investigated if this could be corrected for by changing the heights at which the electric fields change; however, this did not yield satisfactory results. The emission is mainly linearly polarized along the direction of the net force, which for this case is perpendicular to the v × B axis. This results in Stokes Q=I ≃ −1 and U=I ≃ 0. Since the transverse current is much larger than for a typical fair-weather shower due to the strong electric field, the relative contribution of charge FIG. 12. Comparison of stokes parameters with the results of a microscopic calculation for a dual-layered atmospheric electric field, see Table II first column for the structure of the atmospheric electric field.
FIG. 13. Same as Fig. 12 for a triple-layered atmospheric electric field, see Table II second column for the structure of the atmospheric electric field. excess radiation is very small, which is reflected in a much smaller spread in the values of U. The circular polarization of the pulse near the core is small, V=I ≃ 0 because the electric fields in the two layers lie in the same plane.
The observed structure is very different from the field configuration used in Fig. 13. Near the core the signal from the lower layer arrives before the signal emitted from the higher layers since the shower proceeds with the speed of light while the radio signal propagates at a reduced speed due to the finite refractivity of air. This results in a large circular polarization near the core, V=I ≃ 0.5. At a distance of 100 m from the core, due to different traveling distances, the situation is reversed and the signal from the upper layers arrives before that of the bottom layer resulting in a reversed circular polarization, V=I ≃ −0.5. Because of the changing relative importance of the different electricfield layers an intriguing dependence of the linear polarization with distance is observed. Near the core the net polarization is oriented at an angle of 45°with respect to the v × B axis, giving rise to Q=I ≃ 0 and U=I ≃ þ1. At larger distances the radiation from the top layer dominates resulting in a linear polarization normal to v × B (Q=I ≃ −1.0 and U=I ≃ 0.) and vanishing circular polarization, V=I ≃ 0. The intensity shows a strong peak near the core, like is seen for fair-weather events, since for the field configuration of Fig. 13 there is no destructive interference. It is seen that for these rather complicated atmospheric-field configurations the results of semianalytic calculation lie very close to those of the microscopic calculation.

IV. CONCLUSIONS
With a relatively simple parametrization of the structure of the leading plasma cloud in an EAS we are able to recover the main characteristics of the radio emission in an analytic calculation with MGMR3D. The basic structure of the emitted pulse and the complete polarization footprint follow the result of a more complete microscopic calculation. This holds over a wide range of frequencies and some nonobvious structure of the height dependence of the induced currents in this cloud.
On the one hand this gives us a better insight into the physics that is important for understanding radio emission. One finding is that the typical pulse shape, a large peak followed by a long tail of opposite polarity, appears due to the radial dependence of the pancake thickness. In the extreme that the pancake thickness has no radial dependence a bipolar pulse is obtained with comparable strengths for the two polarities. Another interesting point is that one might have expected that the induced transverse currents in the plasma would result in a net dipole charge distribution in the plasma cloud, moving with the speed of light. We have not seen any evidence for such a contribution to the radio pulse. A finite contribution of such a term in the MGMR3D calculation will give rise to a worse agreement with the results of the microscopic calculations.
We are able to rather accurately predict the structure of the radio footprint for rather complicated structures for the height dependence of the induced current using relatively few parameters for the structure of the charge cloud that are kept constant. We see this as a reflection of shower universality. Because of this it is feasible to use MGMR3D in an optimization code to extract the transverse current structure by fitting the results of an MGMR3D calculation to data. To this end the MGMR3D code has been implemented in a Levenberg-Marquardt minimization procedure, which is based on a steepest descent method, to extract the current distribution in the atmosphere during thunderstorms.

APPENDIX A: CALCULATION OF THE RADIATION FIELD
The radiated electric field is derived from the vector potential using A difficulty in evaluating the radiation fields lies in the fact that the retarded distance D ¼ 0 lies in the integration regime. To regularize this one can efficiently use partial integration techniques as shown in Ref. [2]. To evaluate the integration over z 0 in Eq. (2) we follow the same approach as used in Ref. [3] and replace the integral in the z-direction by an integral over λ ¼ h c − h where at the critical height, h c , D ¼ 0. This substitution allows for an easier calculation of derivatives of the vector potential that are necessary to calculate the electric field. When going from the expression for the vector potential to the electric fields the derivatives on the λ integration limits vanish. Note that the retarded time in Eq. (2) depends on λ as well as the other integration variables, as shown in Ref. [2,3].
In calculating the radiation fields we distinguish the charge excess and the transverse current contributions, denoted by the superscripts CX and TC, respectively. The full radiation field is the vectorial sum of the two. The current distributions in the plasma cloud are defined following Eq. (6).

Charge excess
The charge excess in the shower front is given by J Q ðzÞ, as defined by Eq. (15), propagating with the speed of light in the −z-direction thus contributing to the zero and the z component of the vector potential, Eq. (2). Because of the axial symmetry the electric fields are given by where r o ¼ d is the distance from the observer to the shower axis. Substituting the expression for the vector potential, integrating over the spatial extent of the charge cloud the radially polarized radiation field can be written as where r 2 s ¼ ðx 2 s þ y 2 s Þ, the retarded distance D is given by Eq. (5), J Q is the net charge as given by Eq. (15), and the defining Eq. (6) is used to introduce the profile functions fðh; r s Þ [Eq. (8)] and wðr s Þ [Eq. (7)]. The integral is separated into an integration over h, a "ray", a line parallel to the shower axis, of the full shower where special care should be devoted to the point where D ¼ 0, followed by one over the radial extent. In the numerical calculation the results for separate rays, I CX ðtÞ, are stored for a grid of r s and d values, where the latter is the distance of the observer to the single ray.

Transverse current
Along a similar line of reasoning as for the charge-excess radiation, the transverse current radiation field is written as with a similar expression for the component polarized in the y direction. Similar to the notation introduced before, we have J 0x ¼ dJ x =dt r and f 0 ðhÞ ¼ df=dh. The results for rays, 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.

APPENDIX B: PARAMETRIZATIONS
The values of the parameters in the parametrization of the plasma cloud, comoving with the shower, are given in Table III, including the defining equation. All parameters have been determined by fitting the results of a Monte Carlo simulation, and the table specifies the figure where the fit is shown. The only exceptions are the parameters for the longitudinal shower profile, X 0 and γ.

APPENDIX C: PROGRAMMING DETAILS
In the numerical implementation we have exploited the axial symmetry of the current densities as much as possible in order to optimize the running time for the calculation of a single footprint. The integration over the plasma cloud, Eqs. (A3) and (A6), is performed in two steps. In the first step an integration is performed along a single "ray", i.e. Eqs. (A5) and (A7), at a fixed distance from the shower axis, r s , and on a grid of distances from this ray to the observer, r so . Since this is done separately for the charge excess and the transverse current contributions the dependence of the components on the azimuth angles are thus simple sine or cosine functions which will be taken into account at the next stage. These ray integrals are stored on a threefold grid on r s , r so , and time values. The final integral for an observer at a distance, r o , from the core of the shower runs over r s and azimuth angles ϕ s where the distance r so is calculated using straightforward geometry. Thus the radial component (due to charge excess) and the two linearly polarized components [due to transverse currents in the v × B and the v × ðv × BÞ directions] of the radiation field are calculated separately on a grid of observer distances, r o . The calculation of the radiation field for a single antenna can be obtained through an interpolation on the r o -grid and vectorially adding the different contributions. Due to this procedure all angle integrals are almost trivial. The integral for a single ray is relatively expensive since retarded distances have to be calculated with care. At many levels in the calculation an interpolation is necessary to obtain a time-dependent signal for a distance between two grid points. To be able to perform an efficient interpolation it was realized that pulses in subsequent grid points generally have a very similar time structure, however, with a relative time shift. The interpolated pulse can thus be approximated most accurately by averaging the time structures, correcting for the relative time shift and subsequently applying an interpolated time shift. This allowed for using a relatively sparse radial grid.
As an alternative to this procedure we initially used an interpolation of the fourier components of the pulses. For cases where the relative time shifts amounted to half the spacing on the time grid this procedure introduced unrealistically strong high-frequency components for which reason we abolished it.
The full calculation of a radio footprint in MGMR3D takes of the order of 5 CPU seconds, independent of energy and almost independent of the number of antennas, to be compared with approximately a CPU day for CoREAS for the same footprint at E ¼ 10 16 eV and 160 antennas. The time for a CoREAS calculation, for unchanged thinning, increases linearly with energy. The MGMR3D code is available upon request from the authors.