Signal model and event reconstruction for the radio detection of inclined air showers

The detection of inclined air showers (zenith angles θ ≳ 65°) with kilometer-spaced radio-antenna arrays allows measuring cosmic rays at ultra-high energies (E ≲ 1020 eV). Radio and particle detector arrays provide independent measurements of the electromagnetic and muonic shower components of inclined air showers, respectively. Combined, these measurements have a large sensitivity to discriminate between air showers initiated by lighter and heavier cosmic rays. We have developed a precise model of the two-dimensional, highly complex and asymmetric lateral radio-signal distributions of inclined air shower at ground — the “radio-emission footprints”. Our model explicitly describes the dominant geomagnetic emission with a rotationally symmetric lateral distribution function, on top of which additional effects disturb the symmetry. The asymmetries are associated with the interference between the geomagnetic and sub-dominant charge-excess emission as well as with geometrical projection effects, so-called “early-late” effects. Our fully analytic model describes the entire footprint with only two observables: the geometrical distance between the shower impact point at the ground and the shower maximum d max, and the geomagnetic radiation energy E geo. We demonstrate that with this model, the electromagnetic shower energy can be reconstructed by kilometer-spaced antenna arrays with an intrinsic resolution of 5% and a negligible bias.

spaced radio-antenna arrays allows measuring cosmic rays at ultra-high energies (E 10 20 eV). Radio and particle detector arrays provide independent measurements of the electromagnetic and muonic shower components of inclined air showers, respectively. Combined, these measurements have a large sensitivity to discriminate between air showers initiated by lighter and heavier cosmic rays.
We have developed a precise model of the two-dimensional, highly complex and asymmetric lateral radio-signal distributions of inclined air shower at ground -the "radio-emission footprints". Our model explicitly describes the dominant geomagnetic emission with a rotationally symmetric lateral distribution function, on top of which additional effects disturb the symmetry. The asymmetries are associated with the interference between the geomagnetic and sub-dominant charge-excess emission as well as with geometrical projection effects, so-called "early-late" effects. Our fully analytic model describes the entire footprint with only two observables: the geometrical distance between the shower impact point at the ground and the shower maximum d max , and the geomagnetic radiation energy E geo . We demonstrate that with this model, the electromagnetic shower energy can be reconstructed by kilometer-spaced antenna arrays with an intrinsic resolution of 5% and a negligible bias.

Detection of inclined air showers with radio antennas
The detection of inclined air showers with radio antennas has recently been demonstrated by the Pierre Auger Observatory [1]. The radio emission from those showers can illuminate large footprints of several square kilometers, as had been previously predicted via simulations [2], eventually exceeding the footprints measurable by particle detectors [3]. This enables the detection of ultra-high energy cosmic rays (UHECRs) up to the highest energies, e.g., 10 20 eV, with radio antennas, as antennas can be sparsely-spaced such that the instrumented area provides sufficient aperture (> 1000 km 2 sr). The Radio Detector of the upgraded Pierre Auger Observatory [4], which will consist of over 1600 radio-antenna stations covering an area of nearly 3000 km 2 , will routinely detect such events. Other experiments such as the envisioned Giant Radio Array for Neutrino Detection (GRAND) aim to detect inclined air showers from UHECRs as well as from ultra-high-energy neutrinos [5].
Radio antennas measure the electromagnetic radiation emitted by electrons and positrons during their propagation through the atmosphere. While most electrons and positrons are ultimately absorbed in the atmosphere, the radio emission experiences no significant attenuation and can be measured by radio antennas tens to hundreds of kilometers away from the emission region (here approximated as the location of the shower maximum). At such distances, particle detectors mostly measure muons while other particles are absorbed in the atmosphere. Measurements of the showers' muon content suffer from a significant ambiguity between the primary cosmic-ray energy and mass, whereas the strength of the electromagnetic radiation has no significant correlation with the primary mass. Combining both measurements allows us to infer the cosmic-ray mass and test hadronic interaction models [6][7][8][9].
For such studies, an accurate reconstruction of the shower energy from radio measurements is indispensable. For vertical showers, i.e., showers with zenith angles θ < 60°, several signal models of the radio emission at ground to reconstruct the shower energy have been proposed [10,11] and used with experimental data [12,13]. In those models, the radioemission footprints are described based on the macroscopic interpretation of the superposition of two emission mechanisms, the charge-excess (Askaryan) and geomagnetic emission [14,15]. Additionally, the models account for the temporal "Cherenkov" compression of the radio emission, where at a characteristics distance around the shower axis, the emission from the entire longitudinal development arrives almost simultaneously causing an enhancement of the coherent signal [16,17]. This imprints an annulus in the emission pattern, here referred to as Cherenkov ring.
Those models are not applicable for inclined air showers, though, as both the interference between the emission mechanisms and the Cherenkov ring are known to change with the ambient atmospheric conditions in the emission region of the air shower and hence with the zenith angle. Furthermore, for inclined air showers, the radio emission is strongly "projected" onto the ground plane. This projection imprints geometrical early-late effects which disturbs the interference pattern between the both aforementioned mechanisms and becomes significant for zenith angles beyond 60° [18].
Therefore, we present a model dedicated to the description of the radio-emission footprints from inclined air showers. The lateral distribution of the geomagnetic emission is individually described by a 1-dimensional lateral distribution function (LDF) after the emission has been corrected for the aforementioned early-late "asymmetry". With increasing zenith angle, the relative strength of the charge-excess emission decreases and the geomagnetic emission dominates the total signal. Hence, the interference between both emission mechanisms is treated as a further asymmetry to the geomagnetic emission. Using a comprehensive set of CoREAS simulations we found that the shape of the LDF for the geomagnetic emission as well as the asymmetries can be described with a single parameter, d max , the geometrical distance between shower impact position at the ground, in the following also refereed to as "shower core", and the shower maximum. The position of the shower maximum, and hence d max , scales in first order with the zenith angle and in second order with the depth of the shower maximum X max . Lastly, the amplitude of the geomagnetic emission can be described by the geomagnetic radiation energy E geo , i.e., the spatial integral over the energy deposit of the geomagnetic emission at ground. Hence, the entire emission footprint is described by two observables: E geo and d max .
Our model describes the radio emission in terms of energy fluence f [eV m −2 ], i.e., the energy deposit per unit area, in the 30 MHz to 80 MHz band. This frequency band is used by most current-generation large-scale radio detector arrays [19][20][21] and in particular by the Auger Radio Detector. The simulations used here describe the conditions of the Pierre Auger Observatory [22], located near the city of Malargüe, Argentina, in the Southern Hemisphere. This concerns the local magnetic field (strength and orientation), the observation height (altitude), and atmospheric conditions. The adaptability of this model to other conditions, including other frequency bands as many next-generation radio experiments [5,23,24] aim to cover higher frequencies and larger bands, is discussed in Sec. 6.
This article is structured as follows: In section 2, we discuss the simulation sets utilized to derive the radio-emission footprint model and reconstruction algorithm as well as its evaluation. A qualitative description of the radio-emission footprints and their asymmetries is given in section 3. In section 4, we develop the signal model. In section 4.4, the electromagnetic shower energy is reconstructed, and the intrinsic performance is evaluated with simulations of a kilometer-sparse antenna array. Finally, we discuss and conclude in sections 6 and 7.

Simulation and signal processing of the radio-emission footprints of inclined air showers
We use two different sets of air shower simulations, one to develop the model for the radioemission footprints, and one to evaluate the reconstruction of the electromagnetic shower energy with this model. The sets differ mainly in their detector layout, i.e., in the positioning of the observers, and the coverage of the phase space, i.e., the distributions for the showers' energy and arrival direction. For development, we use simulations with an artificial, unrealistically dense detector layout with an antenna grid which is centered around the shower core. The phase space is covered uniformly by discrete, equidistant bins in energy, zenith angle, and azimuth angle. For validation, we use simulations with a realistic, sparse detector layout with showers randomly located within a finite array and a phase space that is sampled continuously in energy and arrival direction. For all air showers, the particle cascades are simulated with CORSIKA [25] and the radio emission is calculated with the CoREAS extension [26]. The detector arrays are placed at an altitude of 1400 m a.s.l. and in a local magnetic field matching the conditions at the site of the Pierre Auger Observatory in Argentina with an inclination of ∼ −36°and a strength of ∼ 0.24 G. If not mentioned otherwise, the simulated atmosphere, i.e., the density profile ρ(h) as function of the altitude h and refractive index at sea level of n = 1 + 3.12 · 10 −4 are used to match the conditions at the Pierre Auger Observatory in October. For the simulation of the particle cascades, an electron multiple-scattering-length factor "STEPFC" of 1 was used. It has been reported previously that lowering this parameter to 0.05 increases the total emitted radiation energy by 11% regardless of the zenith angle and energy [27], but increases the computational effort per shower by a factor of ∼ 4. Hence, we choose to retain a value of 1 with the consequence that the final normalization of our model has to be adjusted for the missing 11% of radiation energy. Below, details to the detector layout for the two simulations sets are given, additional information is summarized in Table 1. For development, we have simulated 4309 showers in which the radio-emission footprint is sampled at 240 observers situated on a flat ground plane such that a star-shaped grid with 8 rays and equidistant antenna spacing is formed in a shower-plane coordinate system perpendicular to the air shower arrival direction (cf. Fig. 1, both panels). Within the showerplane coordinates, the observers are placed depending on the orientation w.r.t. the magnetic field vector to allow for an optimal decomposition of the geomagnetic and charge-excess emission, as it will become clear later.
For evaluation, we have simulated 15970 air showers for which the shower core is randomly distributed within a finite hexagonal array with a spacing of 1500 m. The array resembles that of the Pierre Auger Observatory and extends across nearly 3000 km 2 . For each shower, all observers within a zenith-angle-dependent maximum distance to the shower axis are simulated.
In addition to the 4309 simulations with the star-shaped antenna grid, we simulated three times 216 proton showers with an energy of log 10 (E/eV) = 18.4 and varying atmospheric conditions. These simulations cover the same arrival directions, have 160 observers on a star-shaped grid and a refined particle thinning of thin = 1×10 −6 . The atmospheric conditions match those at the site of the Pierre Auger Observatory in February and June, as well as the US standard atmosphere, as provided within CORSIKA.
In the following, we refer to the atmospheres also with their CORSIKA IDs: US standard: 1, February: 19, June: 23, October: 27. For developing the signal model and reconstruction of air showers, we rely on a model of the atmosphere, i.e., a model for the atmospheric density gradient, from [32,33], which was extended and improved in the context of this work to replicate the atmosphere simulated in CORSIKA/CoREAS.
We use thinning to reduce the considerable computational effort. However, thinning affects the simulation of weak radio signals at large axis distances which have to be treated with caution. This is explained and addressed in more detail in Appendix A.1.
CORSIKA computes a Gaisser-Hillas fit to the energy-deposit table to determine the depth of the air-shower maximum, X max . We found that this fit does not reliably work for air showers with zenith angles beyond 80°. Hence, we perform a simple 2-step χ 2 minimization ourselves to determine the depth of shower maximum. The resulting X max distribution is uncorrelated with the zenith angle [34, Figure 4.1].
The electric field pulses are simulated in the North-South (NS), West-East (WE), and Vertical (V) polarizations. From the time series of each polarization, the energy fluence The emission is simulated for 240 observers (indicated by the gray dots) situated at ground and interpolated in between them using Fourier modes [35]. Left: Largely elongated and strongly asymmetric radio-emission footprint in ground-plane coordinates with the shower incoming from the bottom-right. Right: The same footprint in the commonly used v × B, v × ( v × B) shower-plane coordinates. The white band indicates the contour of 90% the maximum signal. The emission pattern is inconsistent with the typical interference pattern of the both emission mechanisms as it is not symmetric w.r.t. the v × B axis, see details in text.
is calculated by a sum over the squares of the electric field amplitudes in a 100 ns time interval centered around the peak [12, Eq. 1] 2 . The peak is defined as the maximum of the quadratic sum of the Hilbert envelopes from all 3 polarizations. The simulated pulses are band-pass filtered to the 30 MHz to 80 MHz band with an idealized rectangle filter. A frequency resolution of ∼ 100 kHz is ensured by zero-padding the traces sufficiently. For developing our signal model we have to decompose the radio emission at an observer in parts originating from the geomagnetic and charge-excess emission. This is accomplished exploiting the known polarization characteristics of both emission mechanisms and is explained in Appendix A.2. Figure 1 shows the radio-emission footprint at ground from a (simulated) 80°air shower coming from South-East. The color map shows the energy fluence f . The panel on the left shows the footprint on the ground plane which is highly elongated along the shower direction covering a large area with a semi-major axis of ∼ 10 km and a semi-minor axis of ∼ 2 km, and exhibits strong asymmetries. The right panel shows the radio-emission footprint in a shower-plane coordinate system. In this representation, the footprint is more circular. The coordinate system is defined by orthogonal unit vectors pointing into the directions v × B, v × ( v × B), and v, where v is the direction of the primary particle trajectory (i.e., pointing exactly in the opposite direction of the shower axis), and B pointing in the direction of the magnetic field vector which points upwards at the location (latitude) of the Pierre Auger Observatory. This coordinate system is commonly used to display the radio emission from extensive air showers as it highlights the interference between the geomagnetic and charge-excess emission: along the v × B axis the interference is maximal while both emission mechanisms are disentangled along the v × ( v × B) axis; see [36] for a comprehensive review of the emission pattern of the radio emission.

Radio emission from inclined air showers
The footprint in the shower-plane coordinates shows a strong asymmetry (roughly) along the v × B axis (x-axis) which is known to originate from the superposition of geomagnetic and charge-excess emission. However, the footprint is not symmetric w.r.t. the v × B axis as highlighted by the white contour marking 90% of the maximum fluence which is found to be rotated counter-clockwise w.r.t. the v × B axis by 38 • . This inconsistency from the interference pattern of the two emission mechanisms can be explained with the so-called early-late asymmetry. For non-vertical showers, observers at the ground which are below the shower axis will measure the radio emission at an "earlier" stage of the shower development, i.e., they are closer to the point of emission (in this work assumed to be at X max ), than observers above the shower axis. Therefore, the expanding electric field will have a higher intensity and consequently, an early observer will measure a stronger signal than a late observer 3 . Additionally, an early and late observer with equal distances to the shower axis will not have the same off-axis or viewing angle, i.e., not the same angle between the line-of-sight from X max to the observer and the shower axis. An illustration of the differences between an early and late observer is given in Fig. 2. Both effects will introduce an asymmetry in the lateral distribution of the emission which becomes relevant only beyond a zenith angle of 60°a nd increases with the distance to the shower maximum d max and the axis distance of the observer. Correcting for these effects reduces the asymmetry in our simulated radio-emission footprints and restores the known asymmetry pattern from the interference of the geomagnetic and charge-excess emission, as we will see later. A qualitative analysis of the position of the maximum energy fluence in shower-plane coordinates before and after correcting for the early-late asymmetry can be found in Appendix B. Finally, when subtracting the charge-excess emission for the overall emission, we are left with the rotational symmetric geomagnetic emission.
While the overall asymmetry in the radio-emission footprints is dominated by the interference between the two emission mechanisms in the lower half of the zenith-angle range we consider here, the early-late effects constitute the dominant asymmetry for the upper half of the zenith-angle range. Note, that this changes for experiments located at different locations on Earth depending on the strength of the local geomagnetic field.
Besides the asymmetry, the ring-like structure of the temporal Cherenkov compression, i.e., the Cherenkov ring, is visible in the emission pattern. The radius of this ring, i.e., the Cherenkov radius r 0 , can be estimated from the base of a cone with its apex at the shower maximum with an opening angle equal to the Cherenkov angle δ Che (h = h max ). For a point source that is moving with the speed of light β = 1, the radius is where n(h max ) is the refractive index at the shower maximum which is a function of the altitude or height above sea level h.  Figure 2. Illustration of an inclined air shower with an early and late observer at the same early-late corrected axis distance, i.e., with the same off-axis angle. To correct for the early-late effects which cause the asymmetry, we "project" signals measured at ground (indicated with green antenna symbol) along the line of sight from antenna to shower maximum into the shower plane (orange antennas). See text for details.
Recently, an additional "apparent" asymmetry in the radio-emission footprint of very inclined air showers with zenith angle beyond 80 • has been reported [3]. In [37] it is shown that this apparent asymmetry can be explained and resolved by a displacement of the whole radio-emission footprint w.r.t. the Monte-Carlo (MC) shower core. This core displacement is explained by the refraction of the radio emission during propagation in the Earth's atmosphere. Here, we account for it by allowing the core coordinates, i.e., the coordinates of the radio symmetry center, to vary from the MC core. The coordinates are found fitting the lateral signal distribution, cf. Sec. 4.2. The displacement also implies that d max changes if calculated between the shower maximum and the displaced core instead of the MC core. However, the effect to d max is below 1% for all zenith angles and thus ignored in the following.

Model for the radio-emission footprints
To describe the radio-emission footprint from inclined air showers, we first have to remove the early-late asymmetry. In Sec. 4.1, a purely geometrical correction for this asymmetry is formulated and evaluated. With this asymmetry removed, we can determine a parameterization of the shape of the symmetric geomagnetic emission in Sec. 4.2. Our approach to subtract the charge-excess from the geomagnetic emission and describe the interference between the two emission mechanisms, respectively, is described in Sec. 4.3. In Sec. 4.4, the geomagnetic radiation energy E geo is introduced.

Geometrical early-late effects
A description of the early-late effects has already been given in the Sec. 3 and is depicted in Fig. 2. To correct for these effects and eliminate the asymmetry, we "project" the observer positions onto the shower plane perpendicular to the shower axis intersecting with the core and thereby correct their axis distances and energy fluences. This correction assumes the radio emission to expand spherically from a point-like source at the shower maximum with the distance d max to the shower plane, hence the electric field amplitudes scale with the inverse is the distance between an observer and the shower plane, while the unit vector e v points in the direction of the primary particle trajectory. With that, the corrections for the energy fluence f and axis distance r of individual observers are described by the following where the subscript "raw" donates the uncorrected observables. Note that due to the notation of v and B, observers in the positive v × ( v × B) direction are early and have a negative z i coordinate, while observers in the negative v × ( v × B) direction are late with a positive z i coordinate. Fig. 3 (left) shows the radio-emission footprint of the same shower as in Fig. 1 corrected for early-late effects using Eq. (4.2). For this early-late corrected footprint, the symmetry w.r.t. the v × B axis is restored as well as the overall asymmetry is decreased. This allows us in the following to describe the remaining asymmetry solely with the interference of both emission mechanisms.
To evaluate this correction, we have simulated an extra set of 17 showers which have observers on a star-shaped grid in the ground plane (equivalent to the development set) and additional observers on a star-shaped grid situated directly in the shower plane perpendicular to the shower axis. The positions of the observers in the shower plane were chosen such that they correspond to the projected, i.e., early-late correct, positions of the observers in the ground plane. In Fig. 3 (right) the lateral distribution for such a shower with observers, simulated both in the ground and shower planes, is shown. The lateral distribution for the observers in the shower plane (orange circles) has, by definition, no early-late asymmetry imprinted and is much more narrow than the uncorrected distribution for the observers in the ground plane (green squares). The early-late corrected lateral distribution simulated at the ground (blue triangles) shows a good agreement with the distribution directly simulated in the shower plane. In the bottom panel, the ratio between the corrected ground signals and shower-plane signals shows only a slight degradation for large axis distances.
A more quantitative comparison is given in Fig. 4 which presents the ratio between corrected and directly simulated signals across 17 showers with zenith angles ranging from 65°t o 85°as a function of the lateral distance. The axis distance is normalized to the Cherenkov radius r 0 according to Eq. (3.1). As seen in the previous example, the accuracy decreases for larger axis distances. The inset shows a histogram of the presented data. The overall correction is better than to within 5 %.

Lateral distribution of the geomagnetic emission
While the total radio signal exhibits asymmetries, the purely geomagnetic emission is assumed to be rotationally symmetric after any geometrical projection effects have been removed. It can thus be described by a one-dimensional LDF. In [11,21], the LDF for the geomagnetic emission of vertical showers is modeled using a quadratic polynomial in an exponential, i.e., a Gauss curve. This allows one to describe the Cherenkov ring, i.e., the initial rise in energy fluence which is followed by an exponential decay 4 . For more inclined showers, the Cherenkov ring increases in radius, causing a more subtle increase of the emission strength close to the shower axis. In previous iterations of our model, we used a polynomial of the 3rd order in an exponential to account for this more subtle increase. That LDF could describe the region around the Cherenkov ring accurately but decayed too rapidly at larger axis distances, undershooting the simulated signal distribution. Now, to accommodate for this phenomenon The inset shows the same data in log scale. The bottom panel shows the relative deviation between the markers and fitted LDF. The tail of the lateral distribution exhibits a nonphysical flattening due to thinning which is compensated for by setting appropriate uncertainties. and improve the description at larger axis distances, we extend a Gaussian by the addition of a sigmoid. This yields a function f GS with 7 parameters, an amplitude f 0 , and 6 parameters defining the shape of the LDF r fit 0 , σ, p(r), a rel , s, and r 02 : The Gauss-parameters r fit 0 and σ can be interpreted as the position, i.e., radius, and width of the Cherenkov ring. However, it should be noted that r fit 0 does not coincide with the axis distance exhibiting the maximum signal strength, in fact it is slightly larger. This is plausible as the emission pattern is a superposition of the ring-like feature and a decaying exponential function. The exponent of the Gaussian, p(r), is fixed to 2 for axis distances smaller than r fit 0 but can decrease for larger axis distances to accommodate a slower exponential decay. This allows for a better description of the tail of the LDF and was already introduced in [11]. a rel regulates the relative amplitude of the sigmoid term with respect to the Gauss term. The dimensionless parameters s and r 02 define the shape of the sigmoid term. Fig. 5 shows the lateral profile of the geomagnetic emission of an example iron shower. The lateral profile of the geomagnetic emission is obtained by subtracting the charge-excess emission using the concept explained in Appendix A.2, i.e., Eq. (A.2), and after applying the early late correction. While fitting f GS we allow for a shift of the core coordinates to compensate for refractive displacement. Hence, the subtraction of the charge-excess emission and early-late correction are recalculated in each iteration of the fitting procedure. For fitting we use the lmfit python package [38] and a χ 2 minimization.
The lateral distribution of the geomagnetic emission is well-described by f GS . In particular the tail (here at around 1000 m) is more accurately described w.r.t. the previous iteration of our model [39]. However, for even larger axis distances of around 1500 m or more, the LDF does not follow the distribution anymore. This flattening of the simulated distribution is not expected for the coherent radio emission from extensive air showers but is rather the result of thinning, cf. Appendix A.1. Therefore, these signals cannot be trusted. To avoid any bias in the fitting of f GS , we use an uncertainty model for the geomagnetic energy fluence with two terms: a relative contribution of 3%, and a constant value per shower of 10 −4 the maximum geomagnetic fluence f max geo of this shower The latter term ensures (relatively) large uncertainties for weak and potentially thinningaffected signals (cf. the large error bars in that figure). The value of 10 −4 was chosen after a manual inspection of many lateral profiles, the value of 3% was optimized to have a χ 2 /n.d.f.-distribution for all showers with a mean around 1.
While it is no problem to fit an LDF with 7 free parameters (+ 2 core coordinates) to a well-sampled simulated event, in experimental data the signal multiplicity is generally much lower. Furthermore, measured signals are subject to uncertainties and start values for the fit parameters are more uncertain. Hence, it is desirable to reduce the number of parameters, constrain the shape of the LDF to physically reasonable forms, and exploit correlations with shower observables. Here, we investigate the correlation of the shape parameters of Eq. (4.3) (all but f 0 ) with d max . It is worth stressing that this includes an implicit dependency on the zenith angle, atmospheric model, and observation height. First, we fit f GS for all showers with a star-shaped grid. We fix the slope of the sigmoid s = 5 as this ensures that the sigmoid is only dominant within the Cherenkov ring, as desired, and generally simplifies the following procedure: We pick a parameter and parameterize its correlation to d max . Next, we fit all showers again but this time fixing the chosen parameter to its parameterization and inspecting the correlation of the next parameter with d max . We repeat this procedure until all parameters are described by functions of d max . The details of this procedure and all parameterizations are given in Appendix C.1.

Parameterization of the charge-excess strength
So far, we have determined the geomagnetic emission from the simulated pulses using Eq. (A.2) based on the known polarization characteristics of both mechanisms. Thereby, the strength of the charge-excess emission which interferes with the geomagnetic emission in the v × B polarization is estimated from the emission in the v × ( v × B)-polarization. Since the chargeexcess emission and thus the emission in the v × ( v × B)-polarization is relatively weak for inclined air showers, this approach is impractical for the use with measured data as it is difficult to obtain an unbiased estimate of the true emission in the presence of ambient, thermal, Galactic, or anthropogenic noise. Hence, we follow an alternative approach where we define and parameterize a charge-excess fraction to determine the geomagnetic emission. With the following definition for the charge-excess fraction 5 a ce ≡ sin 2 α · f ce /f geo , (4.5) which solely depends on the (dominant) emission in the v × B polarization. The sine of the geomagnetic angle α, the angle between the magnetic field vector and shower axis, accounts for the scaling of the geomagnetic emission with the orientation of the shower to the magnetic field. The cosine of φ, the polar angle between the observer position and the positive v × B axis, accounts for the superposition with the charge-excess emission. Using a parameterization for the charge-excess fraction and Eq. (4.6) rather than Eq. (A.2) has the additional advantage that the parameterization can also be used to subtract the charge-excess emission for pulses close to or on the v × B axis.
In the following, we use CoREAS simulations to derive a parameterization for a ce . First, we extract the charge-excess fraction from the simulated pulses with Eq. (A.2). As mentioned earlier, these equations lose validity for observers close to the v × B axis. Hence, we only consider observers with | cos φ| < 0.9. Furthermore, we select only pulses that are not affected by thinning (cf. Appendix A.1). In Fig. 6 the lateral distribution of the charge-excess fraction for all selected pulses of all showers is shown. The lateral distance is given in terms of the off-axis or viewing angle, and d max is color-coded. The following behavior can be observed: First, the overall strength of the charge-excess emission decreases with increasing distance to the shower maximum, and second, it increases with the lateral distance. The former phenomenon has been studied in simulations for the total energy release between both emission mechanisms in [41], and could be shown in data as well [40]. In contrast, the charge-excess emission increases with the density and hence decreases with the zenith angle, and d max respectively. The scaling of the emission strength of both emission mechanisms with the density at the shower maximum ρ max is discussed in Appendix D. The correlation with the lateral distance has already been reported in Refs. [21,42]. Those observations led to our "ICRC19"-parameterization of the charge-excess fraction [39]: Here, we present a refined version of this parameterization. We substitute the different terms with p ce,i i = 0, 1, 2, as indicated in the formula above. In an iterative procedure p ce,i are optimized fitting Eq. 4.6 with a ce = a ce (p ce,i ) to f GS using the parameterizations established with the procedure presented in the previous section, cf. Appendix C.1. Then, the correlations of p ce,i with d max , ρ max , and r are re-evaluated. The details are given in Appendix C.2. It is worth stressing that ρ max can be determined from d max for a given atmospheric model and zenith angle, and thus does not introduce a new observable/fit-parameter. Finally, we can re-formulate the charge-excess fraction as a function of the r and d max for a given zenith angle, observation height, and atmospheric model: The geomagnetic emission of our example shower, estimated using this parameterization and Eq. (4.6) with the early-late corrected energy fluence in the v × B polarization, is shown in Fig. 7. Compared to the representation in Figs. 1 and 3 (left), the footprint is now fairly rotational symmetric and can be described with the rotational symmetric LDF introduced in Sec. 4.2. Inspecting the footprint closely it becomes apparent that, although the footprint is rotationally symmetric, it is not centered around the coordinate origin which coincides with the MC shower axis. This is due to the refractive displacement of the radio-emission footprint mentioned earlier and described in [37].
It is worth mentioning that significant asymmetries in the lateral distribution of the charge-excess emission were reported in [43] and attributed to shower-to-shower fluctuations. This introduces an irreducible but modest scatter of the charge-excess fraction (see evaluation in the next paragraph). On top of this, an additional dependency on the (azimuthal) arrival direction is apparent in Fig. 17 (first 3 panels), highlighted by the color code, especially for the highest zenith angles (at which the overall relative strength of the charge-excess emission is lowest). Those characteristics are not yet understood and hence not described. They might be related to a so far unexpected dependence of the geomagnetic radiation on the orientation of the geomagnetic field vector which is shown in Sec. 4.4 and further discussed in Sec. 6. However, due to the low relative strength of the charge-excess emission compared to the geomagnetic emission, the remaining scatter does not significantly deteriorate the accuracy as the following evaluation shows.

Reconstruction of the electromagnetic shower energy
So far, we have related the shape of the signal distribution (the symmetric LDF as well as the asymmetry corrections) to d max . What remains is the absolute normalization f 0 . It is easy to see that this parameter correlates with the overall emitted geomagnetic radiation energy E geo , the 2d spatial integral over the f geo footprint at the ground. We can rewrite the LDF to explicitly correlate the signal distribution to E geo with f 0 set to unity. The integral in the denominator has to be solved numerically. The maximum integration distance of 5 r 0 is sufficiently large to evaluate the integral without losing any significant signal. Now we can describe the entire radio-emission footprint with two fit parameters only, E geo and d max (+ two core coordinates). E geo is strongly correlated with the electromagnetic shower energy E em and hence can serve as energy estimator. It should be noted that rather than symmetrizing measured signals f v× B by extracting the geomagnetic emission and applying the early-late correction to compare it to the geomagnetic LDF, it is more practical to apply the asymmetry correction inversely to the LDF to predict the asymmetric signal f v× B which is directly measured and only depends on the shower arrival direction.

Reconstruction of inclined air showers with a sparse antenna array
We use the second set of simulations, with the realistic 1.5 km-spaced antenna array and continuous distributions in arrival direction and energy, to reconstruct E geo and d max with the fully parameterized signal model established in the previous section, and, in a second step, establish the correlation with the electromagnetic shower energy E em . For the definition of E em and how it is determined in simulations see Appendix E.1.
For the following reconstruction, we select only QGSJETII showers with a zenith angle θ > 68 • and at least 5 simulated observers (no requirements on the signal strength of the simulated pulses are imposed). 6210 out of 7972 QGSJETII showers fulfill these requirements. From those 6210 showers we select 6194 showers with a good reconstruction quality. To improve the correlation with E em , we compensate for the second-order scaling of E geo with the geomagnetic angle and local air density at the shower maximum following the logic established in [41], and obtain a corrected geomagnetic radiation energy S geo :  inclined air shower with θ ∼ 75 • . Finally, we can correlate S geo with E em using a power-law: The normalization with ρ has direct implications on the value of S 19 which can be interpreted as the geomagnetic radiation energy for a 10 EeV cosmic ray air shower with an air density at its shower maximum of ρ max = 0.3 g cm −3 . All four parameters S 19 , γ, p 0 , and p 1 are determined in a combined fit of the fitted E geo and ρ max (determined given the fitted d max ) to the true E MC em . Their values are given in Tab. 2. The correlation between S geo and E em (left panel) as well as the achieved reconstruction accuracy (right panels) are shown in Fig. 9. The ratio E em /E MC em is shown once as a function of the true electromagnetic energy E MC em (middle) and once as a function of the true zenith angle (right). The top panels show the full distributions in discrete bins while the bottom panels show the achieved bias (µ) and resolution (σ). The reconstruction accuracy does not depend on the energy with a resolution of better than 5% for all energies. The right panels demonstrate a minor degradation of the energy resolution for the lowest and highest zenith angles. The correlation of E em with the air density, which can be described with the second part of Eq. (5.1), is illustrated in Fig. 10 (left). The y-axis shows y/ y for which has the dependency of the shower energy and geomagnetic angle removed. A significant correlation is visible which is well described by the fitted exponential model, i.e., the second  term in Eq. (5.1). The fitted model agrees well with the one found in [41], although different values for ρ prevent a direct comparison of the fitted parameters. The color code shows the sine of the geomagnetic angle and highlights an unexpected residual correlation which is further discussed in Sec. 6. This residual correlation is partially responsible for the worsening of the energy resolution at larger zenith angle. Fig. 11 shows the ratio E em /E MC em as a function of the true X MC max for each shower (blue dots). The binned mean and standard deviation (error bars) are highlighted by the red markers, the uncertainties on the means are indicated by the error caps. A bias with X MC max is visible: for larger X MC max , the reconstructed electromagnetic energy is underestimated. The overall distributions of X MC max and E em /E MC em are shown as histograms at the top and right axes, respectively. A potential X max -dependent bias in the energy reconstruction is delicate as it could yield a primary-particle dependent bias. However, more than 97% of events are contained within X MC max < 900 g cm −2 for which the bias is below 5% 6 . Furthermore, we did not observe any significant bias in the electromagnetic energy reconstruction between the different primaries. Nonetheless, in a future iteration of this reconstruction this could be improved as discussed in Sec. 6.

Reconstruction of the distance to the shower maximum
In Fig. 10 (right), the reconstructed d max is compared to its true value. The comparison shows an overall good accuracy with no significant bias and a resolution of ∼ 3% which does not significantly depend on the zenith angle. It should be mentioned that the superb reconstruction accuracy of d max achieved here is mainly driven by the fact that we are using the true arrival direction in the reconstruction. In a realistic, experimental setup, where the arrival direction for inclined air showers is only known with a typical accuracy of 0.5° [44] 7 , the reconstruction accuracy of d max will decrease, while we expect that the accuracy on E em is only marginally effected by this. The potential to use d max to reconstruct X max is discussed in Sec. 6.

Reconstruction of air showers generated with a different high-energy hadronic interaction model
We repeated the same evaluation of the air-shower reconstruction with a sparse, realistic antenna array with the Sibyll2.3d-generated air showers. To reconstruct the electromagnetic energy E em , the parameters γ, p 0 , and p 1 are fixed to the values obtained with the QGSJETII showers to allow a direct comparison of the S 19 parameter. S 19 decreased by less than 2% from 3.15 GeV to 3.10 GeV for the Sibyll showers as compared with the QGSJETII showers. It is worth stressing that this change is not due to differences in the prediction of the muonic shower component between both hadronic interaction models, because we reconstruct the electromagnetic energy of the shower. The achieved resolution is very comparable between showers from both interaction models with the exception that the resolution for air showers with zenith angle θ < 70°for Sibyll showed a small degradation in resolution (σ Sibyll Eem (θ < 70°) 7% as compared to σ QGSJET Eem (θ < 70°) 5%). See Fig. 18 in Appendix E.2. The other results, e.g., the reconstruction of d max and the E em reconstruction bias with X max remain practically unchanged.

Discussion
In this section, we elaborate on several features of the here presented signal model and energy reconstruction. In particular, we address the question how the signal model and established parameterizations adapt to different ambient conditions and frequency bands.
The model presented here for the radio emission from inclined extensive air showers in the frequency band of 30 MHz to 80 MHz is tailored to the ambient conditions of the Pierre Auger Observatory. While the general concept and considerations should transform well to other experiments, i.e., other ambient conditions and frequency bands, like GRAND, the explicit parameterizations require revisions. For example, it is known that the Cherenkov ring is more prominent at higher frequencies [37], hence a re-parameterization of the shape of the 1-d LDF seems necessary. Relying on atmospheric models for the parameterization of the charge-excess emission and parts of the LDF parameters reduces the dependency on a particular set of ambient conditions. However, the explicit use of the distance to shower axis r in the parameterization of the charge-excess fraction (4.8) and LDF model (4.3) carries a dependency on the observation altitude.
We normalized our parameterizations to showers arriving perpendicular to the Earth's magnetic field, hence changing orientations of the magnetic field should not affect the model. In [41], the scaling of the (geomagnetic) emission with the strength of the Earth's geomagnetic field was investigated and found to be E geo ∼ B 1.8 . This scaling should apply to our model as well. However, it should be noted, that in Fig. 10 a residual correlation with sin α is apparent which is not yet understood. This correlation becomes stronger with decreasing air density. In [45, Fig. 2] the density scaling of the geomagnetic emission was also investigated but with a stronger geomagnetic field and for the radio emission at higher frequencies. A strong suppression of the geomagnetic emission at lower densities was found, while the behavior at larger air densities is in agreement with the results shown here and in [41]. While we do not see such a suppression here, a common causation related to the magnetic field strength, i.e., the magnitude of the Lorentz force, seems reasonable. In [46], a transition for the geomagnetic radio emission from the regime of time-varying transverse currents to a regime of synchrotron radiation is predicted for air showers developing in low air density in the presence of strong magnetic fields. We do not observe a clear transition in the phase space covered by our simulations (as it is also not expected), however the residual correlation with sin α as well as the suppression of the emission for very low air densities reported in [45] might be explained by this transition. If this explanation proves to be accurate, which needs further investigations, the model presented here requires an adaptation for sites with stronger magnetic fields than present at the site of the Pierre Auger Observatory.
The reconstruction of the electromagnetic energy presented here corresponds to an idealized case and hence the achieved resolution can be considered the intrinsic resolution of the method for air showers reconstructed with a sparse antenna array. In measured data, neither the true arrival direction nor the exact signals arriving at the observers are known perfectly. Ambient and internal noise, an inaccurate detector description (especially of the directional response pattern of the antennas), and other effects will affect the signals reconstructed for each antenna. A detailed study of all those effects is beyond the scope of this paper. A more realistic study of the achievable reconstruction accuracy and reconstruction efficiency with this model has been conducted in [47].
In addition to its application in an event reconstruction, this signal model can also be used to predict the radio emission in the 30 MHz to 80 MHz band from inclined air showers of a given set of energies and arrival directions. This allows studying different aspects of the detection of inclined air showers, for example the effect of the observer spacing on the detection efficiency, when time-and CPU-consuming Monte-Carlo simulations are not available.
Besides the electromagnetic shower energy, the distance of the shower maximum is an observable of great interest as it can be used to determine the slant depth of the shower maximum X max . X max is of special interest as it is commonly used to infer the mass composition of cosmic rays. The distance to the shower maximum d max is reconstructed with a superb accuracy of σ dmax /d max = 3% as shown in Fig. 10 (right). While we intentionally do not probe a realistic scenario, it is worth mentioning (again) that the true arrival direction (zenith angle) is used, with which d max is mostly correlated. Furthermore, small relative changes in d max correspond to large (absolute) changes in X max . Even with a relative d max resolution of 3%, the absolute resolution for the depth of the shower maximum is σ Xmax ≥ 50g cm −2 at d max ≥ 75 km. This leads to the conclusion that the sensitivity of the shape of the lateral signal distribution at ground to X max is rather limited for inclined air showers due to the large distance between shower maximum and detector and the (relatively) small variations in d max induced by variations of X max . While this is unfortunate for obvious reasons, when one wants to estimate the cosmic ray energy it is advantageous as it minimizes the dependency of the LDF to the mass of the cosmic ray primary.

Conclusions
Measuring inclined air showers with radio antennas is of particular interest for two reasons. First, their large footprints allow us to instrument huge areas with sparse antenna arrays, which are necessary to observe the spectrum of cosmic rays at the highest energies. Second, inclined air showers observed in coincidence with radio and particle detectors offer the unique potential to measure the muonic shower component (with the particle detector) and electromagnetic shower component (with the radio detector) independently of each other. The combination of this complementary information yields a strong sensitivity towards the mass of the cosmic ray. For a precise study of the mass composition of UHECR, the energy resolution provided by the radio detector is of critical importance. We present a signal model for the radio emission in the 30 MHz to 80 MHz band from inclined air showers. The model enables accurate reconstruction of the electromagnetic shower energy with a sparse radio-antenna array, an intrinsic resolution of better than 5% and no bias (< 1%) on the primary particle mass. The model relies on an explicit modeling of the dominant, rotationally symmetric geomagnetic emission as well as effects which disturb this symmetric emission and lead to the highly asymmetric pattern we expect from inclined air showers. Those asymmetries are associated with the interference of the charge-excess emission with the geomagnetic emission as well as the imprint of geometrical early-late effects. We exploit correlations between the model parameters and shower observables to minimize the number of free parameters. The final model only relies on two free parameters, the distance between detector and the shower maximum d max and the geomagnetic radiation energy E geo , plus two coordinates for the location of the impact point of the air shower. This allows a reliable fit of the signal distribution and thus efficient event reconstruction.
The presented concept for the signal model is applicable for a variety of radio experiments trying to reconstruct inclined air showers. The described procedure can be used to tune the model parameterizations to match with different ambient conditions as well as different frequency bands relevant for a specific experiment.

Acknowledgments
We would like to thank Alan Coleman for his suggestion to use a sigmoid in addition to a Gauss to describe the lateral profile of the geomagnetic emission as well as for his comments to our manuscript. Also, we would like to thank the anonymous reviewer whose detailed comments helped to greatly improved the quality of this paper. Furthermore, we are thankful to our colleagues involved in radio detection within the Pierre Auger Observatory for very fruitful discussions. We thank our colleague M. Gottowik for his contribution to the simulation library with the PLEIADES cluster at the University of Wuppertal.

A.1 High-frequency emission artifacts from particle thinning
To compute the radio emission from (inclined) air showers with reasonable computational effort, a technique called thinning is used [48]. This implies that particles produced in a single interaction and below a certain energy threshold are removed from the simulation except for one randomly selected particle. This particle is assigned a weight to describe the ensemble of particles it "replaces" such that energy conservation is preserved. The probability for a particle to be selected is proportional to its energy. This dramatically reduces the number of particles to be simulated while correctly reproducing showers on average. Random particle fluctuations and thus shower-to-shower fluctuations are affected. However, if the energy threshold E th = thin E 0 and the maximum weight a particle can be assigned w max = thin E 0 /GeV which both depend on the parameter thin , are chosen wisely [48], the effect is tolerable. For the simulation of the radio emission, thinning introduces another problem. A particle with a large weight, which represents many particles, emits a radio wave with an amplitude scaling with its weight (with no phase difference) while the actual ensemble of particles emits radio waves with phase differences. In other words, particles which are described by one particle with a corresponding weight emit perfectly coherent emission. This effectively introduces artificial additional power. For small lateral observer distances, this power is well below the actual coherent radio emission. However, for increasing lateral distances or when considering higher frequencies, i.e., with decreasing coherence, this artificial signal starts to significantly impact the simulated power and subsequently, the affected pulses need to be rejected.
In the left panel of Fig. 12, the spectra of two pulses are presented. The observer of one pulse is closer to the shower axis (top) and the other one further away (bottom). For both pulses, the spectra of the v × B-and v × ( v × B)-polarizations are shown, representing the signals of the geomagnetic and charge-excess emission contributions, respectively, as both observers are situated along the v × ( v × B) axis. The band of interest from 30 MHz to 80 MHz is highlighted. Both spectra show the same feature: A smooth exponential decay of the amplitude followed by a noisy plateau. While the first is expected for coherent emission, the latter is not and thus interpreted to be caused by thinning. While the pulse of the closer observer is not (or in the case of the v × ( v × B)-polarization only slightly) affected by the noise floor in the band of interest, the pulses of the observer further away from the shower axis show a significant disruption in both polarizations and thus have to be rejected from further analysis. To quantitatively examine whether a pulse is contaminated or not, we fit a first-order polynomial to the logarithmic spectrum in the frequency ν range between 30 MHz to 80 MHz, i.e., |A(ν)| = 10 mν ·ν+b (A.1) with a slope parameter m ν and a constant b. The slope parameter m ν as a function of the lateral distance for an example shower is shown in Fig. 12 (right). While the spectrum is almost flat (m ν = 0) on and around the Cherenkov ring, it is falling more steeply with increasing lateral distance as expected. Around 750 m a kink is visible. The lateral distance of the observer whose pulse in the v × B polarization has the steepest slope after which the disruption in the considered band becomes considerable, defined as r min , is identified. To be conservative, we select pulses only from observers with a lateral distance smaller than r thin min = 0.85 r min as clean, pulses of observers with larger lateral distances are considered affected by thinning artifacts. For the example event, the dashed line indicates this criterion. The considered maximum lateral distance per shower scales in first order with the zenith angle and just slightly with the energy. This is expected since energy-dependent weight limitation was used [48] to simulate the air showers. For highly inclined showers, observers with lateral distances of over 2 km are still considered. With this selection, the number of considered observers decreases from 240 to around 160 per simulated shower. This selection is solely used for the parameterization of the charge-excess fraction in Sec. 4.3 given the fact that it is otherwise difficult to independently identify affected observers and mitigate their effect on the parameterization. For fitting the lateral distribution of the geomagnetic emission (cf. Sec. 4.2), we consider all observers -even the ones with pulses affected by thinning artifactsbut assign an appropriate uncertainty to all signals, effectively reducing the impact of weak signals, to avoid any bias from affected observers.

A.2 Decomposition of the radio signal
For an accurate description of the radio-emission footprints it is useful to decompose the emission into the geomagnetic and charge-excess contributions which is possible due to their polarizations characteristics. The geomagnetic emission is polarized in the negative v × Bdirection while the charge-excess emission is polarized radially inwards [36]. First, the electric field traces simulated in the [ e NS , e WE , e V ] coordinate system are rotated into the [ e v× B , e v×( v× B) , e v ] coordinate system. This allows us to calculate the energy fluence for each of these polarizations f v× B , f v×( v× B) , and f v , while f v is almost zero since the electric field of the radio emission is oscillating perpendicular to v. Then we decompose the signal into one part originating from the geomagnetic f geo and another part originating from the charge-excess f ce effects, making use of the known polarization characteristics, i.e., (derived from [11]) while the term 1/ sin 2 φ diverges, hence the ansatz loses validity. It should be noted that this ansatz is not affected by the early-late asymmetry as this asymmetry does not affect the polarization of the emission. However, the disentangled signals f geo and f ce need to be early-late corrected to show the expected symmetry. A different ansatz that overcomes this problem but comes with other disadvantages is discussed in [34, Appendix B.1] and only mentioned here for completeness. It has to be mentioned that the equations in (A.2) assume that both emissions arrive simultaneously at an observer, i.e., without any phase shift. Such a phase shift would give rise to a circularly polarized component in the incoming electric field which indeed has been seen in experimental data [49], i.e., there is a time delay between the pulses originating for the charge-excess and geomagnetic emission. To quantify the fraction of circular polarization in the radio pulses we calculate the Stokes parameters I, Q, U , V following the procedure detailed in reference [49]. Since the relative strength of the charge-excess emission decreases with the zenith angle (cf. Sec. 4.3), the fraction of circularly polarized signal is small for most showers in our set. The determined time delay, following [49], is within ∆t < 1 ns for most observers and thus the above equations are applicable for the radio emission in the 30 MHz to 80 MHz band. To ensure that this holds, we only use showers with a geomagnetic angle α > 20 • , i.e., the angle between the shower axis and the Earth's magnetic field vector, to develop the model.

B Effect of the early-late asymmetry on the emission pattern
As explained in Sec. 3, in addition to the asymmetry due to the interference between geomagnetic and charge-excess emission, there is also an early-late asymmetry present in the radio-emission footprints from inclined air showers. The latter disturbs the well-known emission pattern produced by the interference from which the maximum emission is expected at the v × B axis. In Fig. 13 the position of the maximum emission is shown in polar coordinate for uncorrected emission patterns (right panel) and early-late corrected emissions (left panel). In both panels, the outer axis gives the rotation from the v × B axis in degree, the inner axis gives the lateral axis distance normalized to r 0 and the color code shows the d max of the respective showers. In the figure, only showers coming from the South are shown as for those the shower axis projected on the ground is in the same plane as the v × ( v × B) axis and perpendicular to the v × B axis. As it can be seen, with increased d max the maximum is rotated towards the v × ( v × B) axis, i.e., towards the incoming direction of the air showers. This indicates that the early-late asymmetry is increasingly dominant. If the early-late effects are corrected, the maximum is found around the v × B axis which is again consistent with the interference pattern from geomagnetic and charge-excess emission only. It is also apparent that for the corrected emission pattern the maximum for showers with a zenith angle of 85 • (= d max ∼ 150 km) is found rotated from the v × B axis. However, keep in mind that at those inclinations the charge-excess emission is vanishing and hence no clear maximum at the v × B axis is expected.

C Parameterizations for the signal model
In the next section C.1, the six parameters of the lateral distribution function f GS (4.3) describing the shape of the geomagnetic emission are correlated with d max . In section C.2 the parameterization of the charge-excess fraction a ICRC19 ce (4.7) is optimized to obtain the new expression a ce (4.8).

C.1 Parameterizations of the shape of the geomagnetic emission
First, we investigate how the radius of the Gaussian r fit 0 relates to d max . In Fig. 14, the opening angle of a cone originating at d max , δ fit Che (r fit 0 ) = tan(r fit 0 /d max ), calculated from the fitted radius r fit 0 , is shown (top panel) as a function of d max (dots). The prediction for the Cherenkov angle δ Che according to Eq. (3.1) is shown for comparison (lines). Both, δ fit Che and the theoretical prediction are shown for four different simulated atmospheres. The atmospheres at the location of the Pierre Auger Observatory for February (summer), July (winter), and October correspond to the maximum, minimum, and yearly average for the refractivity at ground level, respectively [32]. The bottom panel shows the relative deviation between fitted and predicted angles. The comparison shows an overall remarkable agreement for larger zenith angles and different atmospheres. For lower zenith angles, a systematic deviation can be found. However, it is possible to use r 0 determined as a function of d max according to Eq. (3.1) instead of fitting it without losing significant accuracy. We carefully checked that the remaining free parameters sufficiently compensate for the deviations introduced when using the predicted value of r 0 . In the following, we refer to r 0 as the Cherenkov radius.
Next, we study the correlation of σ in Eq. We normalize the function with the term "d max −5 000 m" to decrease the statistical fluctuations in the fitted parameters. However, this restricts the parameterization to values of d max > 5 000 m. While d max < 5 000 m is very unlikely for hadron-initiated air showers with zenith angles θ ≥ 60 • as it would require depths of X max > 1200 g cm −2 , for neutral particles, in particular neutrinos, d max < 5 km is not difficult to imagine. However, if we assume that the radio emission is still detectable at a 4 • off axis angle, the maximum axis distance for d max = 5 km is ∼ 350 m, which is too small to detect such showers in more than one or two antennas with a kilometer-spaced detector. The uncertainties of the fitted data are statistical, estimated from the χ 2 -minimization of the LDF fit. They can not explain the deviation of single data points from the parameterization. It can not be excluded that those points represent an alternative minimum. However, the global minimum can be easily identified with the parameterization by the vast majority of the data. To obtain the optimal values for the parameterizations in Eq. (C.1) we employed again a χ 2 -minimizations, this time using the iminuit python package [50]. The same procedure is now applied consecutively to the parameters p(r) (resp. b), a rel , and r 02 , in this order. Their distributions are shown in Fig. 15 and their parameterizations are given by Eqs. (C.2) -(C.4): In the distributions for a rel and r 02 , an additional trend, not described by the parameterizations, is significant. Within one zenith angle bin, a steep increase of the corresponding parameter from deep to shallow showers is apparent. The matter is further discussed in the last paragraph of this section, for now we choose to only describe the correlation of all parameters with d max .
We also verified the fit results for different atmospheres. With the prediction of r 0 depending on the atmospheric profile, the parameterization of the LDF f GS explicitly uses information of the atmosphere. The other parameters, however, are assumed to be universal, i.e., do not depend (significantly) on the simulated atmosphere. In Fig. 16, the correlation of the parameters with d max for the different simulated atmospheres is shown. Although the atmosphere influences the correlations of the parameters with d max , the variation is minimal and the October atmosphere used for the parameterization indeed describes the mean reasonably well. The effect on the geomagnetic radiation energy was found to be below 1%. However, it should be kept in mind, that this conclusion was obtained with star-shaped simulation which for this purpose has limited informative value.
The observed X max -dependent bias in the energy reconstruction can be resolved, to a large degree, by describing the secondary correlation of X max with the LDF parameters a rel or r 02 , cf. Fig. 15. The secondary correlation can be explained by the ambiguity of d max for different zenith angles and X max values. This ambiguity, although not completely resolved, can be reduced using the density at the shower maximum ρ max as observable. An elegant solution to resolve the ambiguity is the introduction of d 750 the distance between shower core and a fixed slant depth of 750 g cm −2 , in the parameterizations of a rel or r 02 . However, when doing so, we found an implausible kink in the distribution of the fitted d max distribution with simulations from the validation set (cf. Fig. 10, right), and for this reason, the fact that we do not observe a significant primary-particle dependent bias and for the sake of simplicity we decided to not include such a term in the parameterizations. However, if needed in the future, our model can be improved by a more thorough study of these secondary correlations.

C.2 Refined parameterization of the charge-excess fraction
To refine the parameterization of the charge-excess fraction a ICRC19 ce (4.7), we optimize the different terms p ce,i in the parameterization by fitting the distribution of symmetrized signals, i.e., signals for which the charge-excess emission has been subtracted using the parameterization, to the rotational symmetric, fully-parameterized LDF from Sec. 4.2. First, we optimize the term describing the scaling of the charge-excess fraction with the air density p ce,2 . While optimizing Eq. 4.6, i.e., f geo (a ce (p ce,2 )), to minimize the difference to the constrained f GS , we find the value for p ce,2 for which the estimated geomagnetic emission is the most symmetric. Only the normalization f 0 of f GS is varied as well. Fig. 17 (Top-left) shows the correlation of p ce,2 obtained for all showers with ρ max . The purple curve shows our new description given by p ce,2 = ρ max 0.428 kg m −3

3.32
− 0.0057. (C.5) The functional form is rather ad-hoc but describes the data better than the exponential function used at the ICRC19 [39] which is shown by the orange curve. Also, this new function can become negative at small ρ max , and thus implausible, but does so for lower values ρ max < 0.09 kg m −3 than the exponential model used in the ICRC19-parameterization which becomes negative for ρ max 0.15 kg m −3 . This allows us to extend the parameterization to zenith angles of 85°and beyond. Similarly, the "exponential correction" p ce,1 term and the "off-axis angle" p ce,0 term are substituted and refined (in this order). For p ce,1 the mean is used. In Fig. 17 (Top-right) the total exponential term is shown using r = r 0 per shower and compared to the previous value. For the off-axis angle term, instead of a constant factor a linear model with a slope depending on d max is used now, cf. Fig. 17 (Bottom-left).

D Density scaling of the geomagnetic and charge-excess emission
While the air-density scaling of the overall emitted radiation energy has been investigated earlier [41] and found to coincide with the scaling of the dominant geomagnetic emission, the correlation of the charge-excess emission with the air density has not previously been studied that thoroughly. We extract the charge-excess emission by Eq. (A.2) from star-shaped simulations and perform a 2d-spatial integration over the interpolated footprint to estimate its radiation energy. We interpolate the 2d-footprint via Fourier decomposition [35]. We find that the charge-excess emission decreases in absolute strength (and not only relative strength) with increasing zenith angle, i.e., decreasing density at the shower maximum, see Fig. 17 (bottom-right). A similar correlation is reported in [45, Fig. 2]. An explanation for this phenomenon might be that in a denser atmosphere, i.e., for more vertical showers, more electrons are ionized from the ambient atmosphere hence the negative charge excess is stronger. However, this simple explanation needs verification. The scaling of the geomagnetic emission with the density of the shower maximum has already been shown in Fig. 10 (left). It can be explained with the following picture: The emission strength depends on the mean free path length with which the electromagnetic particles traverse the atmosphere. With a larger mean free path, equivalent to traversing a less dense atmosphere negatively and positively charged particles can drift further apart before interacting, creating a stronger transverse current and thus resulting in a stronger geomagnetic emission. For a given slant depth, the density at the shower maximum ρ max , is smaller for larger zenith angles and larger d max , respectively.

E.1 Deriving the true electromagnetic shower energy from CORSIKA simulations
The strength of the radio emission is strongly correlated with the energy of the electromagnetic shower cascade, i.e., the electromagnetic shower energy. It should be stressed that this is slightly different from the fluorescence light seen by optical telescopes which better correlates with the total calorimetric energy (for which other particles like muons have a non-negligible contribution). We compute the electromagnetic energy as the sum over the longitudinal energy deposit E i for gamma rays, electrons, and positrons (ionization and cut) as provided in the CORSIKA DATnnnnnn.long files, i.e., It is worth noting that this includes the energy deposit in the ground plane (which is accounted for in the two last rows of this table with the SLANT option). In inclined air showers, no clipping effects of the radio emission occur, because the showers can evolve fully before the ground is reached.

E.2 Reconstruction of the electromagnetic shower energy for showers generated with Sibyll2.3d
Fig . 18 shows the E em reconstruction for showers generated with the Sibyll2.3d high-energy interaction model. From the 6199 showers with a zenith angle greater than 68°and at least 5 simulated observers, 6185 showers were reconstructed with good quality. The results are very comparable to the ones archived with the QGSJETII-04 showers: The S 19 parameter decreased slightly by less than 2% and the resolution at lower zenith angles worsens slightly. It is worth stressing that the simulations used to develop this model were solely generated with QGSJETII-04, hence a small decrease in reconstruction quality for Sibyll-generated showers is not surprising. Nevertheless, this result underlines the fact that the radio emission has little dependence on the underlying hadronic interaction model as long as one normalizes quantities to the electromagnetic energy in the air shower.