An analytic description of the radio emission of air showers based on its emission mechanisms

Ultra-high energy cosmic rays can be measured through the detection of radio-frequency radiation from air showers. The radio-frequency emission originates from deflections of the air-shower particles in the geomagnetic field and from a time-varying negative charge excess in the shower front. The distribution of the radio signal on the ground contains information on crucial cosmic-ray properties, such as energy and mass. A long standing challenge is to access this information experimentally with a sparse grid of antennas. We present a new analytic model of the radio signal distribution that depends only on the definition of the shower axis and on the parameters energy and distance to the emission region. The distance to the emission region has a direct relation to the cosmic ray’s mass. This new analytic model describes the different polarizations of the radiation and therefore allows the use of independently measured signals in different polarization, thereby doubling the amount of information that is available in current radio arrays, compared to what has been used thus far. We show with the use of CoREAS Monte Carlo simulation that fitting the measurements with our model does not result in significant contributions in both systematic bias and in resolution for the extracted parameters energy and distance to emission region, when compared to the expected experimental measurement uncertainties. This parametrization also enables fast simulation of radio signal patterns for cosmic rays, without the need to simulate the air shower.


Introduction
Ultra-high energy cosmic rays (UHECRs) impinging onto the atmosphere induce huge cascades of secondary particles. Established techniques for their detection are the measurement of the particles of the air shower that reach the ground, the observation of the isotropic fluorescence light emitted by molecules that have been excited by the shower particles [1,2] or by non-imaging air-Cherenkov telescopes that measure the incoherent Cherenkov light produced by the shower particles [3]. Important observables for most analyses of high-energy cosmic rays are their energy and the atmospheric depth of the shower maximum X max , which is an estimator of their mass. In particular the accuracy, i.e., the systematic uncertainty, of the energy measurement is a crucial aspect. The determination of the cosmic-ray energy from stand-alone particle detectors needs to rely on Monte Carlo simulation, where the hadronic interactions have large uncertainties. So far, the best accuracy is achieved with the fluorescence technique, but this is only possible at sites with good atmospheric conditions. Furthermore, precise quantification of the scattering and absorption of fluorescence light under changing atmospheric conditions requires extensive atmospheric monitoring efforts [1].
Another independent method for the detection of cosmic rays is the detection of broadband radiofrequency emission from air showers [4,5]. The radio technique combines a duty cycle close to 100% with an accurate and precise measurement of the cosmic-ray energy [3,6,7,8], as well as a good sensitivity to the mass of the primary cosmic-ray [9]. In particular, the energy measurement is well-compatible with, and may even outperform, the fluorescence technique in terms of achievable accuracy [10,11]. This is mostly due to the transparency of the atmosphere to radio waves and the corresponding insensitivity to changing environmental conditions, and because the radio-frequency emission can be calculated theoretically via first principles from the air-shower development [12,13].
The radio emission from air showers is due to the acceleration and creation of charged particles within the air shower [14] and is described by classical electrodynamics. In practice, particles other than electrons and positrons do not contribute significantly to the radio emission due to their smaller charge-to-mass ratio [4]. From a macroscopic point of view, radio emission is attributed to two main emission mechanisms: The geomagnetic and charge-excess emission processes. In the dominant geomagnetic emission process, electrons and positrons are deflected in the geomagnetic field in opposite directions due to the Lorentz force, resulting in a transverse current. The strength of the emission scales with sin α, where α is the angle between the particle movement (shower axis) v and the geomagnetic field B. In the charge-excess emission process, a time-varying negative charge-excess in the shower front leads to a longitudinal current which is mostly due to the knock out of electrons from air molecules.
The spatial distribution of the energy fluence, i.e., the energy per unit area of the radio electric-field pulse, holds information on relevant air shower parameters such as the energy and the atmospheric depth of the shower maximum X max [15]. The amount of energy emitted in the form of radio emission by the air shower -referred to as the radiation energy -is given by the spatial integral over the energy-fluence. The radiation energy is directly related to the electromagnetic shower energy E em and allows for a precise measurement with a theoretical energy resolution of only 3% [12]. Thus, the radiation energy serves as a universal estimator of the cosmic-ray energy and is already exploited by the Pierre Auger Collaboration to measure cosmic-ray energies [7,8].
The shape of the spatial signal distribution is primarily determined by the distance D Xmax from the observer to the emission region. The emission region can be approximated by the position of the shower maximum X max [12]. The distance D Xmax depends primarily on the zenith angle θ of the air shower and scales approximately with D Xmax ∝ 1/ cos θ, with a second order dependence on the value of X max for the typical physical range of X max [12]. The dependence is visualized in Fig. 1 left. The usage of D Xmax has the advantage that a universal description of the radio signal distribution can be given that does not depend on the specific altitude of the experiment.
A long-standing challenge to access the energy and X max information experimentally with a sparse grid of antennas is an analytic modeling of the radio signal distribution and will be addressed in this article. In [16], an empirical parametrization for the spatial radio signal distribution is introduced based on morphological arguments, which gives an adequate description of the data measured by LOFAR and the radio array of the Pierre Auger Observatory (AERA) and was already successfully exploited to measure cosmic-ray energies [7,8]. However, explaining the behavior and value of the parameters of this parametrization is not straightforward, as most parameters depend on various shower features. With the knowledge gained over the past years (e.g. [12,16,17]), we formulate an analytic description of the spatial signal distribution directly based on its physical emission processes whose parameters directly depend on the air-shower parameters energy, incoming direction and X max . In addition, we explicitly use the polarization of the radio signal which effectively doubles the available information of each antenna station. This is achieved by the following approach: We model the spatial signal distribution on the ground originating from the geomagnetic and the chargeexcess emission separately. Then, the two signal-strength distributions are both radially symmetric around the shower axis [12]. We note that for inclined air showers an additional asymmetry due to the projection of the signal distribution on the ground arises. This imposes no principle problem for our approach but requires an additional correction of the projection effect first. Hence, we restrict our analysis to air showers with zenith angles smaller than 60 • where the projection effect is still negligible. Then, the asymmetric two-dimensional radio signal distribution is modeled naturally by the interference of the two emission mechanisms. This is because the two emission mechanisms exhibit distinct polarization signatures. The geomagnetic emission is polarized in the direction of the Lorentz force v × B acting on the shower particles. The charge-excess emission, in contrast, is polarized radially towards the shower axis. The parametrization presented here, also enables the fast simulation of expected signals in a radio detector array. Starting with the species and energy of an incoming cosmic ray and a choice of X max and direction of the cosmic ray, the antenna signals can be predicted for any antenna position relative to the shower axis from simple geometric considerations.
In the following, we first present the Monte Carlo data set that we used to develop an analytic description of the geomagnetic and charge-excess function. Then, we present the geomagnetic and charge-excess functions separately and exploit the correlations of the parameters of the functions with the air-shower parameters. Finally, we combine the two functions to model the two-dimensional radio signal distribution. Throughout this work we follow the maxim of practical usability of this function, i.e., we demand a precise description of the data with a sufficiently small number of parameters so that it can be applied to current radio air-shower detectors. Following this maxim, we also offer a reference implementation in python that is available on github [18].

Monte Carlo data set and decomposition of radio signal into geomagnetic and charge-excess contributions
We use the CoREAS program [19], the radio extension of the CORSIKA code [20], for the simulation of the radio-frequency emission from air showers. In CoREAS, each shower particle is tracked and the radiation resulting from its movement is calculated from first principles using classical electrodynamics [14]. The radio emission originates only from the movement of electrons and positrons as the contribution from e.g. muons is negligible due to their higher mass-to-charge ratio. This allows for a precise calculation of the radio emission as the development of electromagnetic showers is well understood.
Recently, many tests have been performed to investigate the accuracy of the radio predictions. On the experimental side, the LOFAR radio cosmic-ray detector with hundreds of antennas with small spacings allows for precise tests of the theoretical predictions. No significant deviation of the CoREAS calculation from experimental data was observed [21,9]. In addition, a detailed comparison of the CoREAS code with the independent ZHAireS [22] air-shower simulation code was carried out and showed no significant difference in the shape of the radio signal distribution [13]. Hence, we can use the CoREAS code to develop a precise description of the radio signal distribution and study the dependence on air-shower parameters.
For this analysis we use a set of 300 air showers simulated with CoREAS 7.5602 with QGSJetII-0.4 [23] and UrQMD [24] as hadronic interaction models. The geomagnetic field is set to an inclination of -35.9 • and a strength of 0.24 Gauss which corresponds to the value at the Pierre Auger Observatory. We note that this choice does not reduce the general applicability of our results. The scaling of the radio signal with the geomagnetic field is well understood [12] and our results can be rescaled to different geomagnetic field configurations. The amplitude of the geomagnetic component scales almost linearly with the magnetic field strength whereas the charge-excess component is unaffected by the magnetic field. In [12] the proper formulas are given to rescale the simulated energy fluences, and we implemented the rescaling in our reference implementation. The thinning level is set to 1 × 10 −6 with optimal weight limitation and the lower energy thresholds for electrons/positrons and photons are set to 250 keV. We use the monthly average atmospheric profile for October at the Pierre Auger site that is available in CORSIKA and corresponds to the yearly average at that site.
A fraction of 50% of the air showers have an iron primary and 50% have a proton primary. The cosmic-ray energy is distributed between 10 17 eV to 10 19 eV, uniformly in the logarithm of the energy. The zenith angle θ is distributed uniformly in cos θ from 0 • to 60 • and the azimuth angle is chosen randomly. For each air shower, we calculate the radio emission for two observations planes, one at 1564 m a.s.l. -corresponding to the altitude of the radio array of the Pierre Auger Observatory (AERA) -and another one at sea level -corresponding to the altitude of the LOFAR cosmic-ray radio detector. A suitable coordinate system is in the shower plane (the electric field is always polarized perpendicular to its direction of propagation v) where one axis is aligned to the v × B direction (the polarization of the geomagnetic component) and the other axis to the v × ( v × B) direction. In each observation plane, the observer positions are positioned in a star pattern in this v × B coordinate system projected on the ground plane (see Fig.1 right). This choice of antenna positions allows for an effective sampling of the two-dimensional radio signal distribution and decomposition of the emission into its geomagnetic and charge-excess contributions (cf. Fig. 1 right and [12,16] for more information about the choice of observer positions). The radio pulses are filtered to be limited to the 30-80 MHz band, which corresponds to the bandwidth of most current cosmic-ray radio detectors [8].
As CoREAS is a microscopic Monte Carlo code, no emission mechanism is explicitly modeled. Therefore, the contribution of the geomagnetic and charge-excess emission processes to the simulated electric field can not be differentiated. However, we can exploit the different polarization signatures of the two emission mechanisms to decompose the signal into its geomagnetic and charge-excess contribution [12].
In Fig. 1 right, the distribution of the energy fluence is shown in the v × B-v × ( v × B) coordinate system for a typical air shower. At observer positions on the v × ( v × B) axis, the polarizations of the signals from the geomagnetic and charge-excess processes are orthogonal. Here, the v × B component of the electric field E v× B originates only from geomagnetic emission, whereas the v × ( v × B) component of the electric field E v×( v× B) originates only from charge-excess emission. Hence, we calculate the geomagnetic and charge-excess energy fluences f geo and f ce from the respective electric-field components: where ε 0 is the vacuum permittivity, c is the speed of light in vacuum and ∆t is the sampling interval of the electric field E( r, t), which depends on the position r, φ (here in polar coordinates) and time t. The positions for φ = 90 • correspond to positions along the positive v × ( v × B) axis (cf. Fig. 1 right). The shape of the spatial distribution of the energy fluence depends on the distance D Xmax from the observer to the shower maximum. We observed three different categories of shapes corresponding to air showers (A) that hit ground before emitting most radiation energy; (B) that hit ground shortly after emitting all radiation energy; and (C) that have large distances between the ground and the air-shower development.
In Figs

Shapes of the signal distribution
The shape of the spatial radio signal distribution depends on the distance between the observer and the emission region. We call this quantity D Xmax and measure it in grammage along the shower axis. This behavior is illustrated in the sketch of Fig. 2 showing the three typical shapes. A characteristic number is D Xmax = 430 g/cm 2 at which the radio emission of the air shower is almost completed (99% of the radiation energy has already been emitted) [12]. This means that if a detector measures an air shower at this D Xmax , the shower development will just have completed when hitting the observer. For this or smaller distances to the shower maximum, the distribution of the energy fluence is peaked and narrow around the shower axis (example A). Thereby, the geomagnetic signal distribution is always narrower than the charge-excess signal distribution. E.g., in Fig. 3, the geomagnetic distribution shows a sharp peak at the shower axis, whereas the charge-excess distribution already flattens at the shower axis and shows a Gaussian like shape.
The second example (B) is for an intermediate distance D Xmax = 572 g/cm 2 , where the shower development is already completed but the observer is not yet far away from the emission region. Now, the geomagnetic signal distribution also flattens at the shower axis and its shape is Gaussian like. In the chargeexcess case, it starts to become visible that the energy fluence drops to zero at the shower axis. This is an expected behavior as the polarization flips, i.e., changes by 180 • , at the shower axis. Only if the energy fluence drops towards zero at the shower axis, do we get a continuous transition. We note that also for smaller distances to X max the charge-excess energy fluence becomes zero at the shower axis, but on such small scales that it is not visible in the finite sampling of our simulations [12].
The third example (C) is for a distance D Xmax = 1046 g/cm 2 , far away from the emission region. In particular, the change of the shape of the signal distribution between the second and third example is due to free propagation of the electromagnetic waves and not because additional radio emission is created by the air shower. For large distances to X max the emission is peaked in a Cherenkov cone, which originates from the refractive index of air being larger than unity. The opening angle of the cone depends on the air pressure at the emission region, i.e., on the height of the emission. The peaking structure of the Cherenkov cone is smeared because the emission occurs in an extended lateral and longitudinal region along the shower axis.
Independent of D Xmax , we observe that the charge-excess component shows more fluctuations than the geomagnetic component, e.g., in examples B and C, the right wing of the signal distribution shows a slightly higher maximum amplitude than the left wing. For other showers in our data set, both wings have the A B C Figure 2: Illustration of how the shape of the distribution of the energy fluence changes with distance to the shower maximum. We note that we only show vertical air showers here for illustration purposes. Throughout the analysis we used air showers with all kinds of incident directions, and the distance to the shower maximum depends strongly on the zenith angle, e.g. example A is a typical shape of a vertical shower whereas example C is a typical example of an inclined shower. same maximum amplitude or the left wing has a higher maximum amplitude than the right wing. Similarly, the upward fluctuation near the shower axis of example A (cf. Fig. 3 top right) vanishes for other air showers or appears at a different position. Hence, this behavior is likely to be attributed to shower-to-shower fluctuations. Accordingly, our goal in this article is to model the underlying smooth signal distribution and not to model single fluctuations, although we recognize that once the underlying signal distribution is well modeled these fluctuations may provide interesting additional information on the shower development of an individual event.

Discussion of measuring D Xmax in grammage vs. geometric distance
In this article we measure the distance from the observer to the shower maximum in grammage (g/cm 2 ) and not in units of the geometrical distance (km). This choice is not necessarily obvious because the width of the function should be a function of the geometric distance to the shower maximum if the observer is far away from the emission region. This is because at large distances the shower development has finished and the radio emission propagates freely through the atmosphere. So one can think of a cone that gets wider the further away the observer is. However, using the geometric distance to the shower maximum comes with the following disadvantages: The transition point, where the shower is fully developed when hitting the surface, has only the same D Xmax for all zenith angles if we measure the distance in g/cm 2 . Hence, we can only describe the transition of the function between the different shapes correctly if we measure D Xmax in g/cm 2 . As a consequence of measuring D Xmax in g/cm 2 , the D Xmax dependence is not completely universal but depends on the model of the atmosphere used in the CORSIKA simulation. Nevertheless, our model will still describe the data/simulations of different atmospheres but the 'D Xmax fit parameter' has a slight offset to the true D Xmax in the order of 10 g/cm 2 -20 g/cm 2 (cf. Sec. 6.2 for more details).
The shower development itself depends on D Xmax as a consequence of the air pressure profile: the larger the air density the shorter the geometrical distance in which the shower develops and thereby the smaller the size of the emission region, which in turn influences the spatial distribution of the radio frequency emission. Due to the change in refractivity with air density over the path from the emission region to the detection plane, the radius of the Cherenkov ring also depends on D Xmax , in addition to depending on the geometric distance of the emission region to the detector. Hence, the spatial signal distribution at the detector depends both on D Xmax and on the geometric distance to X max . Although these dependencies are not the same, they are rather similar and tend to be degenerate in a fit, except when a huge number of measurements at a large variety of positions are available.
As the goal of this article is to describe the radio signal distribution over the complete D Xmax range, we chose to measure D Xmax in g/cm 2 . For a model of the signal distribution dedicated to horizontal air showers, i.e., only for high zenith angles where the observers are far away from the shower development (case C of Fig. 2), one could use the geometric distance to the shower maximum or a mix of D Xmax and geometric distance. The latter choice would probably also remove the second order dependence of the signal shape on the observation height (cf. Fig. 8).

Signal distribution of the geomagnetic emission
The strength of the geomagnetic emission is circularly symmetric around the shower axis and, thus, only a function of the perpendicular distance to the shower axis r. In the v × B coordinate system, r is given by The energy fluence of the geomagnetic emission can be parameterized as The parameter R geo can be interpreted as the radius of the Cherenkov ring, and the parameter σ geo describes the width of the function. For R geo > 0, the function can be interpreted as signal from a smeared Cherenkov ring that contributes from both sides of the shower axis and thereby fills up the central area in a natural way. The function p(r) is a small correction to an exponent of 2 and will be discussed below. For p(r) = 2, the two-dimensional integral over the function, which gives the radiation energy, can be calculated analytically. The constants N R− and N R+ are chosen such that the parameter E geo is the geomagnetic radiation energy for p(r) = 2. A visualization of the function is shown in Fig. 6a. Negative values of R geo describe the situation when the air shower has not yet emitted all radiation energy when hitting the observer. Then, the signal distribution is strongly peaked around the shower axis and is described by the falling flanks of a Gaussian function (cf. Fig. 3). Positive values of R geo describe the distribution of the energy fluence after the shower has emitted all its radiation energy (which is roughly at D Xmax ≈ 430 g/cm 2 [12]). Then, the function is the sum of two Gaussian functions centered at R geo and -R geo . If the radius R geo becomes larger then the width σ geo , the function becomes peaked at the Cherenkov ring.

Determination of optimal parameters
In this section, the parameters of the geomagnetic function that describe the simulations best are determined in a χ 2 minimization. For each simulated air shower, we do not use all data points in the fit but exclude data points with energy fluences smaller than 10 −4 of the maximum energy fluence where the simulation shows large fluctuations and may be less reliable. The challenge in this multi-parameter fit is to find a procedure such that the global minimum is found correctly for all signal shapes. We therefore employ the following procedure: The variation of the exponent is fixed to p(r) = 2 in the first iteration. As the geomagnetic radiation energy can also be calculated via a numerical integration of the data points, we also fix E geo in the fit to the result of the numerical integration and determine the optimal parameters σ and R geo in a χ 2 minimization. In the fit, the data points with large energy fluence are given a larger weight than those with small energy fluence to have the central part of the function described well. This is done by assigning the same absolute uncertainty to all data points. The fit result can be seen as a dashed blue curve in Figs. 3 -5. For larger distances of the observer to X max (examples B and C where R geo is positive), we observe that the central part of the signal distribution is described well, but the signal fall-off at large distances is slightly overestimated (cf. left panels of Fig. 4 and 5). This can be modeled by a modification of the exponent of the exponential function of Eq. (3) of the following form: The functional form of p(r) is visualized for typical values of r cut and b in Fig. 6b. For positive values of b, p(r) becomes smaller than 2 for distances larger than r cut . Hence, the signal fall-off at large distances weakens. We determine the optimal parameters of r cut and b again in a χ 2 minimization, where we fix all other parameters (E geo , σ geo and R geo ). We give all data points the same relative uncertainty to increase the weight of the data points with small energy fluence. With this modification the simulated signal distribution can be described well at all distances r. The geomagnetic function with the modification of the exponent is shown as a solid orange curve in the left panels of Fig. 4 and 5. For small distances to X max , where the signal distribution is modeled with negative values of the parameter R geo , a variation of the exponent p(r) is not necessarily required. However, we observe that the parameters R geo and σ geo are not very well constrained. Larger negative values of R geo can be compensated by larger values of σ geo such that the χ 2 function has no clear minimum. While this imposes no problem for the reconstruction of the radiation energy, the correlation between R geo and D Xmax is disturbed. We found that allowing for a variation of the exponent p(r) solves this problem. As the parameters r cut and b are correlated with σ geo and R geo , a separate fit as in the R geo > 0 case will not work. Instead, we first determine E geo , σ geo and R geo with p(r) = 2 as described above with the additional constraint of R geo > −200 m. Then, we fit all five parameters of the function simultaneously with the start parameters of E geo , σ geo and R geo set to the values of the previous fit result and without any constrains on R geo . In the combined fit we assign the same relative uncertainty to all data points to increase the weight of the data points with small signal strength.

Dependence of fit parameters on the distance to the shower maximum
The goal of this section is to reduce the number of fit parameters by exploiting correlations with the distance to the shower maximum. In the end, the geomagnetic function should depend only on the radiation energy, which is the integral of the function, and D Xmax which determines its shape. In the following, we first parametrize the variation of the exponent p(r, r cut , b) (Eq. (4)) as it is only a small correction to the shape of the function. Then, we repeat the fit with p(r, r cut , b) = p(r, D Xmax ) fixed to its D Xmax parametrization (note that D Xmax is a known quantity in this simulation study). This will result in a more stable determination of the optimal parameters of σ geo and R geo and less fluctuations in their correlation with D Xmax .
Parametrization of p(r). In the left panels of Fig. 7, the correlations of r cut and b with D Xmax are shown for all air showers in our data set. The correlation shows a different behavior for functions with R geo > 0 and with R geo < 0. In addition, some fits converged at large negative values of R geo , which result in b parameters close to zero. Therefore, we ignore all data points with R geo < −250 m in the parametrization of b(D Xmax ). 1 The correlation of r cut and b with D Xmax are both parameterized with spline functions. Spline functions have the advantage of being capable of describing arbitrary relations analytically with a small set of parameters. The technical details of using spline functions are discussed in Appendix B. In how much detail a spline function describes the data depends on its number of parameters. This smoothness is controlled by an external parameter during the determination of the optimal spline function. We adjust the smoothing  condition manually such that the main trend of the correlation is followed but the function does not follow a single fluctuation. In particular, we adjusted the smoothing condition such that the parametrization has a smooth transition between the R geo < 0 and R geo > 0 cases. The resulting spline function is shown as green curve in Fig. 7. The parameters of the function are tabulated in our reference implementation [18].
Parametrization of σ geo and R geo . In the next step, the geomagnetic function (Eq. (3)) is fitted again to our data set with p(r) fixed to its D Xmax parametrization. Then, the correlation of σ geo and R geo with D Xmax can be studied, which is presented in the left panels of Fig. 8. We observe that the correlation is slightly different for observers at different altitudes. Hence, we perform a separate parametrization for an observation altitude of 1564 m a s l (the height of the Pierre Auger Observatory) and at sea level (the height of the LOFAR detector). Again, the flexibility of spline functions allows us to parameterize the correlations in a continuous and smooth way. Their parameters are presented in [18]. In case of fully developed air showers (D Xmax 430 g/cm 2 ), both σ geo and R geo show a smooth, nearly linear increase with D Xmax . The parameter R geo increases faster than σ geo . At around 600 g/cm 2 , R geo becomes larger than σ geo resulting in a visible Cherenkov ring. For smaller D Xmax , the dependence is more complex and difficult to interpret due to the interplay between σ geo and R geo . In particular, it is difficult to model the transition from a Gaussian shaped to a narrowly peaked signal distribution, i.e., the transition from fully developed showers to showers that are still developing when hitting the observer. Although the individual dependencies σ geo (D Xmax ) and R geo (D Xmax ) are not monotonous, their combination leads to a geomagnetic function f geo (D Xmax ) that is smooth in D Xmax . This is, with increasing D Xmax , f geo shows a smooth transition from a narrowly peaked distribution, via a Gaussian shaped distribution, to a broad distribution with a visible Cherenkov ring. Thus, it fulfills the primary objective of this article: An analytic description of the radio signal distribution whose shape is determined by one variable only, the distance to X max . As this behavior is not directly obvious from Fig. 8, we provide a video of the development of the geomagnetic function with D Xmax as supplemental material [25].
We now managed to parameterize the geomagnetic function in terms of only two air-shower parameters: E geo and D Xmax . The only issue that still needs our attention is that E geo does not correspond directly to the radiation energy E geo because of the variation of the exponent p(r). This is because Eq. (3) can be integrated analytically only for p(r) = 2 and we normalized the function only for the p(r) = 2 case (see Appendix A). However, we can integrate Eq. (3) numerically for each value of D Xmax and parameterize the deviation between E geo and E geo as a function of D Xmax . Again, we use splines to parameterize this relation and present the parameters in [18]. The final function is presented later in the text in Eq. (7).

Signal distribution of the charge-excess emission
The strength of the charge-excess emission is circular symmetric around the shower axis and can be described with a modification of the Gamma distribution with k ≥ 0. The variation of the exponent p(r) has the same functional form as in the geomagnetic case (cf. Eq. (4)). For k = 0 (and p(r) = 2), the function is a Gaussian function with mean zero. For k > 0 the function has the property to be zero at r = 0. The interplay between the rising part from r k and the falling part from the exponential function models the Cherenkov ring. The distance where the function becomes maximal is given by R ce = σ ce √ k/ √ k + 1. The constant N ce is chosen such that the two-dimensional integral over f ce is E ce for p(r) = 2. Hence for p(r) = 2, E ce equals the radiation energy of the charge-excess emission E ce . A visualization of f ce is shown in Fig. 6c.
For small distances to X max , the signal distribution is maximal and peaked at the shower axis, which is described with k = 0 (cf. Fig. 3). We note that for k = 0 the energy fluence does not become zero at the shower axis. Here, the electric-field vector changes its sign and the energy fluence should become zero. However, as discussed above, the drop towards zero at the shower axis occurs at such small scales that it is not detectable in any experiment as the typical size of an antenna is larger than the distance at which the energy fluence becomes zero. In particular, the sampling of our simulations is not fine enough to see this effect at small distances to X max . Also in practice, not modeling this effect in our function does not have any effect on the radiation energy (the integral over the function) nor on its dependence on D Xmax . For larger distances to X max , it becomes visible that the energy fluence goes to zero at the shower axis. Hence, k becomes larger than zero to model the observed behavior (cf. Fig. 4 and Fig. 5). We find that a modification of the exponent p(r) leads to better results at large distances to the shower axis for all distances to X max .

Determination of optimal parameters
To obtain the optimal fit parameters, we follow the same procedure as for the geomagnetic case. We first determine the parameters σ ce and k in a χ 2 minimization where the radiation energy is fixed to the result of a numerical integration of the data points, p(r) is fixed to 2, and all data points are given the same absolute uncertainty to increase the influence of the data points with high energy fluence. The resulting charge-excess functions are shown as blue dashed curves in Figs. 3 -5. Then, in a separate fit, the optimal parameters r cut and b of the variation of the exponent are determined. In this fit, E ce , σ ce and k are fixed to the previous fit results and the same relative uncertainties are given to all data points.

Dependence of fit parameters on the distance to the shower maximum
As in the geomagnetic case, we first parameterize the relation of b and r cut with D Xmax , which is shown in the right panels of Fig. 7. We observe a different behavior for k > 0 and k = 0. For k = 0, the parameter r cut is always zero and the dependence of b on D Xmax can be described with a straight line. For k > 0, the relation between b and D Xmax is also described by a straight line. However this time, r cut is not zero but shows a correlation with b that can be described by a second order polynomial. From this correlation we can calculate the dependence of r cut on D Xmax The parameters of the parametrization are presented in [18] and in Appendix C.
In the next step, the charge-excess function is fitted again to our data set with p(r) fixed to its D Xmax parametrization. Then, the correlation of σ ce and k with D Xmax can be studied, which is presented in Fig. 8. In the case of σ ce , we again observe that the correlation is slightly different for observers at different altitudes. Hence, the correlation is parameterized separately for 1564 m a s l and sea level. We model the correlation with B-spline functions and present their parameters in [18].
For the parameter k we do not observe any difference between different observation altitudes. Before the air shower has emitted all its radiation energy (at D Xmax ≈ 430 g/cm 2 ), k is zero. For larger distances D Xmax , k increases monotonously with D Xmax . The relation can be described well with a logistic function of the form The optimal parameters are listed in [18] and in Appendix C. We provide a video of the development of the charge-excess function with D Xmax as supplemental material [25].

Extrapolation to larger zenith angles
In this analysis we considered only air showers with zenith angles up to 60 • . This is because an additional asymmetry becomes relevant for larger zenith angles due to the projection onto the ground -which we did not take into account in this work. The reason for this asymmetry is that at different observer positions the radio signal traversed different amounts of atmosphere until it reaches the ground. Observer positions 'below' the shower axis see smaller distances than positions 'above' the shower axis resulting in a left-right asymmetry which becomes relevant above 60 • zenith angle.
However, our results indicate that the parameters of the function increase monotonically with increasing D Xmax and we assured ourselves that our parametrizations follow the observed trend to at least D Xmax = 2000 g/cm 2 (the corresponding zenith angles can be read off from Fig. 1 left). In particular for the observation altitude of 1564 m a s l , the extrapolation to larger D Xmax are similar to the 0 m a s l simulations that have ∼400 g/cm 2 larger D Xmax values at 60 • zenith angle (cf. Fig. 1 left and Fig. 8). Hence, our results can likely also be used at larger zenith angles (up to D Xmax = 2000 g/cm 2 ) if the additional asymmetry due to the projection effect is taken into account. However, to make this model usable for horizontal air showers in general, the parametrization of the parameters of our function with D Xmax should be extended to larger D Xmax values using new CoREAS simulations.

Combination to two-dimensional function
With the results of the last two sections, we are able to describe the geomagnetic and charge-excess energy fluence distributions with only three parameters: The radiation energies of the two emission processes and the distance to the shower maximum. We can reduce the number of parameters further by noting that E ce can be expressed as a function of E geo using the result of [12]: The relative charge-excess strength is a function of the air density at the shower maximum. With a model of the atmosphere, i.e., a description of the density as a function of height, the density at X max can be calculated from the distance to X max . Then, we can express both functions as a function of the complete radiation energy E rad = E geo + E ce and D Xmax .
Then, the geomagnetic function reads where c geo , a, R geo , σ geo and p(r) are functions of only D Xmax and c geo (D Xmax ) is the parametrization of the ratio between E geo and E geo . The parameter a is the relative charge-excess strength and sin α is the angle between the shower axis and the geomagnetic field. Similarly, the charge-excess function reads where c ce , a, σ ce , k and p(r) are functions of only D Xmax and c ce (D Xmax ) is the parametrization of the ratio between E ce and E ce .
To obtain the total energy fluence at any position we need to combine the radially symmetric geomagnetic and charge-excess functions and take the interference between the two components into account. For the electric field, the simple relation holds at any position. We can write down this relation explicitly for the two components of E: Fig. 1 right and describes the position relative to the shower axis. From this relation and Eqs. (1) + (2), the interference in units of the energy fluence can be calculated [12]: This calculation assumed that the geomagnetic and charge-excess component are in phase which introduces a small overestimation of f v× B of 1% as studied in [12]. In Sec. 6.1 we do not see that the resolution of the radiation energy and D Xmax is negatively impacted by this approximation and 1% is anyway much smaller than the typical uncertainty on the energy fluence in an experiment of e.g. 5% in case of the radio array of the Pierre Auger Observatory [8]. A phase difference can straightforwardly be introduced at the expense of an extra parameter by multiplying the E ce component in Eq. (9) by a factor cos(ζ) with 0 ≤ ζ ≤ π the relative phase difference between the geomagnetic and charge excess component. Such a phase parameter could be included in a fit to experimental data.
Using these relations, the geomagnetic and charge-excess energy fluences can be combined to the total observed energy fluence at any position. In Figs. 9-11, the total energy fluence of our three example air showers at different distances to X max are presented with the optimal fit result of the combined two-dimensional function. We note that the function is completely defined by the two parameters E rad and D Xmax , i.e., only these two parameters are optimized, and that the minimization is very stable once the function is parametrized to depend only on these two parameters. The energy fluence is shown along the v × between the two components. Along the v × B axis, the two components are polarized into the opposite direction for negative distances and are polarized into the same direction for positive distances. Hence, we observe a destructive interference on one side of the shower axis and a constructive interference on the other side of the shower axis. This demonstrates that the observed asymmetry along the v × B axis is modeled well by the interference between geomagnetic and charge-excess emission. To better picture the evolution of the total observed energy fluence with D Xmax , we again provide a corresponding video as supplemental material [26].

Precision of analytic description
In this section we quantify how well the function describes the simulated signal distribution. However, there is no unique way to do this. E.g., comparing the relative difference of energy fluences at each position will mostly highlight differences at large distances to the shower axis where the absolute difference is small. And calculating absolute differences of energy fluences normalized to the maximum of the function -the method that was used in [16] -will mostly highlight (dis)agreement at the maximum energy fluence. Although describing the maxima of the distribution with high precision seems like the most important thing, it is not for most analyses. The physical quantities that are extracted from the distribution of the energy fluence is the radiation energy and D Xmax , and both of these quantities depend little on a precise modeling of E rad = 9.7e+02 MeV D Xmax = 1044 g/cm 2 geo+ce LDF geo ce Figure 11: Same as Fig. 9 but for the air shower of Fig. 5.
the maximum amplitude. For determining the shape, it is more important to model the falling part of the distribution correctly. And for the radiation energy, the parts of the distribution where the product of energy fluence f times the distance to the shower axis r is largest are most important, because of the larger covered area. Therefore, we judge the precision of the analytic description by how well the radiation energy and D Xmax is extracted from the distribution of the energy fluence. This is done by using the parametrization that depends only on E rad and D Xmax , i.e., Eqs. (12)- (14) together with (7) and (8). This function is fitted to the two-dimensional distribution of the energy fluence and the fitted values of the radiation energy and distance to X max are compared with the true MC values. The outcome of this study is presented in Fig. 12. The radiation energy can be determined with a resolution of 4% and D Xmax with a resolution of 13 g/cm 2 . We note that with the knowledge of the zenith angle of the air shower D Xmax can be converted to X max , which is an estimator of the cosmic-ray mass. As these values are much smaller than the typical experimental uncertainties on E rad and D Xmax that originate mostly from a finite sampling of the energy fluence and uncertainties in the measurement of the energy fluence itself, our analytic description is sufficiently good and does not limit the experimental resolution of sparse radio arrays. Even for detectors with a high station density, such as LOFAR, our model can serve as a fast alternative to the computationally intensive template matching technique [9] without dominating the X max resolution.
We developed this function with the prime goal of usability at radio arrays and shared our work with the Pierre Auger and the LOFAR collaboration right from the beginning. Therefore, this function has already been tested for sparse radio arrays and has even been successfully applied to data from the Pierre Auger Observatory and the LOFAR cosmic-ray detector [27,28] to determine cosmic-ray properties.

Usage at different observation altitudes and atmospheric conditions
Some of the parameters of our function do not only depend on D Xmax but also show a slight dependence on the observation altitude as presented in Fig. 8. Therefore, we presented separate parameterizations of R geo , σ geo and σ ce for an observer at sea level and at 1564 m a s l . As the shapes of the parameterizations for the two observation altitudes are very similar, a pragmatic way to use this function for intermediate observation heights is to linearly interpolate between the two parameterizations. Inspecting Fig. 8 also allows to estimate an upper limit on the uncertainty for a specific observation height. To get a better assessment of the uncertainties, the fit results of this adjusted model can be compared to a few CoREAS simulations produced for the new observation altitude.
This analysis was performed for a specific profile of the atmosphere that corresponds to the yearly average at the Pierre Auger Observatory. Different atmospheric conditions result in a change of D Xmax , e.g., the D Xmax values increase by 10 g/cm 2 for small D Xmax and up to 20 g/cm 2 for large D Xmax values if the atmospheric profile is changed to the US standard atmospheric model. As a consequence the D Xmax parameter of our function (called D fit Xmax in the following) does not correspond to the true value D true anymore. We note that the function still describes the data well and that the inferred radiation energy is accurate, it is just that the fit parameter D fit Xmax will differ from the true D true Xmax value. Hence, depending on the accuracy that is required, one can easily calculate the relation D mymodel Xmax (D thismodel Xmax ) [29] using only the two atmospheric models. To achieve a better accuracy, new CoREAS simulations for a particular atmospheric model need to be produced. Then, our model, as it is, can be fit to the new simulations to determine the relation D new MC Xmax (D fit Xmax ). With this prescription, different experiments do not need to go through the exercise of reparametrizing the 'low level' fit parameters (σ geo , R geo , σ ce , k) with D Xmax , and the number of required simulations will be less than the ones used in this work.

Conclusion
Ultra-high energy cosmic rays can be measured through the detection of radio-frequency radiation from air showers. The radio emission originates from deflections of the air-shower particles in the geomagnetic field and from a time-varying negative charge excess in the shower front. The distribution of the radio signal on the ground contains information on crucial cosmic-ray properties, such as energy and mass. The strength of the radio emission scales with the cosmic-ray energy and the shape of the spatial signal distribution depends on the distance to the emission region. A long standing challenge to access this information experimentally with a sparse grid of antennas is a corresponding analytic description.
We have presented a new analytic model of the radio signal distribution that models the spatial distribution of the energy fluence originating from the two emission mechanism separately. The observed two-dimensional asymmetry in the signal distribution is modeled by the interference between the two emission mechanisms. Thereby, we explicitly take into account the polarization of the radio signal by separately describing the v × B and v × ( v × B) component of the energy fluence. Hence, the available information at a single antenna station is doubled which allows for a more precise determination of the signal distribution compared to previous models for the same number of detector stations.
One parameter of our model is the radiation energy, which directly relates to the electromagnetic shower energy. The contribution to the resolution of the radiation energy of our model is 4% which translates to an uncertainty of the cosmic-ray energy of only 2% due to a quadratic scaling between the two quantities. This is negligible to practical sampling uncertainties of real radio detector arrays. Hence, our model is particularly useful to precisely determine the cosmic-ray energy from radio air-shower measurements.
The remaining parameters of our function are correlated to the distance to the shower maximum. We have shown that this model can determine the distance to shower maximum to a precision of 13 g/cm 2 . The experimental sampling limitation of a real radio array will likely dominate the final uncertainty also here. In addition, there will be some deterioration when converting D Xmax into X max , especially at larger zenith angle, due to the zenith angle resolution of the shower and atmospheric conditions uncertainties. However, the fitting procedure presented here is unlikely to be the limiting factor in the final X max resolution for practical radio arrays. Hence, we can formulate our model to depend only on the radiation energy, the distance to X max and the core position, which was always at the coordinate origin in this simulation study. Thus, our model provides direct access to the main air-shower parameters energy and X max and manages to use the smallest possible number of parameters. Furthermore, our model provides an absolute prediction of the energy fluence at any position for a given air-shower energy and X max . For many studies our model can replace computational-extensive simulation studies, e.g., for estimating the sensitivity of a detector. We provide a reference implementation of our model in Python [18].
For the charge-excess function the dependence of r cut on D Xmax is described by The parameters that define the various spline functions are presented in the reference implementation [18] and in the supplemental material [30].