Refractive displacement of the radio-emission footprint of inclined air showers simulated with CoREAS

The footprint of radio emission from extensive air showers is known to exhibit asymmetries due to the superposition of geomagnetic and charge-excess radiation. For inclined air showers a geometric early-late effect disturbs the signal distribution further. Correcting CoREAS simulations for these asymmetries reveals an additional disturbance in the signal distribution of highly inclined showers in atmospheres with a realistic refractive index profile. This additional apparent asymmetry in fact arises from a systematic displacement of the radio-emission footprint with respect to the Monte-Carlo shower impact point on the ground. We find a displacement of $\sim 1500\,\text{m}$ in the ground plane for showers with a zenith angle of $85^\circ$, illustrating that the effect is relevant in practical applications. A model describing this displacement by refraction in the atmosphere based on Snell's law yields good agreement with our observations from CoREAS simulations. We thus conclude that the displacement is caused by refraction in the atmosphere.


Introduction
Radio detection of inclined air showers is a powerful technique for the detection of ultra-high energy cosmic e-mail: felix.schlueter@kit.edu † e-mail: gottowik@uni-wuppertal.de rays. It has been demonstrated that the detection of these particles is possible with sparse antenna arrays [1]. The radio emission, which originates from the electromagnetic component of an air shower, experiences almost no attenuation while propagating through the atmosphere and hence provides an accurate and precise energy estimator [2]. Combined with a measurement of the muonic shower component, e.g. by a particle detector, measurements of inclined air showers also provide a high mass-composition sensitivity [3], which will be in particular a target of the large-scale radio detector of the upcoming upgrade of the Pierre Auger Observatory [4]. Therefore, inclined air showers have a high relevance and a detailed understanding of the signal distribution of their radio emission is crucial to accurately reconstruct the cosmic-ray properties of interest.
The radio signal distribution is affected by a strong asymmetry arising from the superposition of the geomagnetic and charge-excess emission caused by their individual polarization patterns [5]. For inclined showers with a zenith angle larger than 60 • an additional early-late asymmetry becomes relevant [6].
In [7] a previously unknown apparent asymmetry in the radio-emission footprint of inclined CoREAS simulations was observed which was presumed to be related to refraction of the radio emission in the atmosphere. In this paper we explain and resolve this apparent asymmetry with a systematic offset of the radio-emission footprint with respect to the Monte-Carlo (MC) air shower impact point, or Monte-Carlo "core". We express this offset as a displacement of the radio-emission center from the Monte-Carlo core. We develop a method to determine the radio core without implying detailed knowledge of the lateral distribution function (LDF) of the radio-emission footprint. Furthermore, we present a model successfully describing the core displacement by the refraction of electromagnetic waves propagating through a refractive atmosphere based on Snell's law.
In [8] the influence of the refractive index on the radio-emission footprint and thus the reconstructed depth of the shower maximum X max was studied. An important correlation of the reconstructed X max with the refractivity at the shower maximum was found. Focusing on vertical showers, effects on signal propagation were, however, neglected.
The effect presented here is crucial in order to gain a more detailed understanding of the signal distribution of radio emission from inclined air showers. Implications arise for the modeling and reconstruction of air showers as well as the interpretation of the reconstructed geometries. Our investigation mainly refers to the frequency band of the radio emission from 30 to 80 MHz. This frequency band is used by most current-generation large scale radio detector arrays [9][10][11] as well as the radio detector of the upgrade Pierre Auger Observatory [4]. However, many next generation radio experiments [12,13] aim to cover higher frequencies and a larger band, e.g. 50 to 200 MHz for the GRAND experiment [12]. In Sec. 3.4 we will therefore address the comparability of our results with this frequency band.
This article is structured as follows: In the following section we show that the additional apparent asymmetry in the lateral distribution of inclined showers found in [7] can be explained with a displaced radio core. In Sec. 3 we present a method to determine the core of the lateral distribution of the radio emission. Using a set of CoREAS showers, we establish a systematic displacement of the radio core with respect to the MC core. Furthermore, we investigate differences in the Cherenkov radius for geomagnetic and charge-excess emission. In Sec. 4 we present a model based on Snell's law that successfully describes the displacement, and we discuss the treatment of the propagation of radio emission in CoREAS. Finally, we draw our conclusions in Sec. 5.

Apparent asymmetry in the signal distribution of air-shower radio emission
For vertical air showers, several LDF models exist which take the interference between geomagnetic and chargeexcess emission into account [14,15]. It is assumed that the geomagnetic and charge-excess radiation, independently, have a rotationally symmetric emission pattern around the shower axis. However, [2] reports relative deviations in the energy fluence of the subdominant charge-excess contribution of 20 %. In inclined air showers the radio emission above the shower axis propagates considerably longer through the atmosphere than below axis. Consequently below the shower Fig. 1 Lateral distribution of the radio emission of a 85 • air shower along the positive and negative v × (v × B) axes with respect to the MC core position. The energy fluence is expected to be symmetric on both axes. This is fulfilled for a constant refractive index n (blue lines), independent of its exact value. If the refractive index changes with height (orange lines), the LDF is not symmetric with respect to the MC shower axis. On closer look, a displacement of the symmetry axis rather than an asymmetry is observed. Figure updated from [7]. axis, i.e., early in the shower, observers at ground measure a higher intensity than observers late in the shower and thus an early-late asymmetry is imprinted in the radio-emission footprint. A geometrical correction consisting of a "projection into the shower plane" assuming a point source of the radio emission located at the shower maximum X max can resolve this asymmetry within 2 % [16].
In addition to these nowadays well-known asymmetries, a further apparent asymmetry was observed in [7]. To further investigate this finding, we simulated one air shower with a zenith angle of 85 • arriving from South for four different atmospheric refractivity profiles using the CoREAS code [17]. In Fig. 1 the lateral distribution of the radio signal along the v × (v × B) axis is shown in terms of the energy fluence f with the unit eV/m 2 (v: incoming shower direction, B: magnetic field vector (B points upwards with an inclination of ∼36 • ) ≡ positive v × (v × B) axis is early in shower). Along this axis the geomagnetic and charge-excess contributions are decoupled by their polarisation [2] and after a correction for geometrical early-late effects no asymmetry is expected. If the refractive index is set to a constant value of n ≡ 1 or n ≡ 1.000 03 (approximately the value of n(h max ) at the shower maximum for an air shower with a zenith angle of 85 • ) throughout the atmosphere, the signal distribution in the shower plane is symmetric along the positive and negative v × (v × B) axes with respect to the MC shower axis. The exact value of the refractive index is not important for symmetry, but changes the shape of the LDF. With a changing refractive index following the density gradient in the atmosphere, an apparent asymmetry is observed, the LDF is not symmetric anymore with respect to the MC shower axis. This is enhanced when doubling the refractivity n − 1 throughout the atmosphere. It seems that for these simulations the symmetry axis is displaced from the MC shower axis in the direction of the positive v × (v × B) axis which is in early direction of the shower. In the next section we will show that this apparent asymmetry in fact originates from a displacement of the radio-emission footprint with respect to the MC core.

Displacement of the radio-emission footprints
We now analyse the apparent asymmetry introduced by the refractive index profile of the atmosphere in detail using a set of 4308 inclined showers simulated with CORSIKA [18] and its CoREAS extension [17]. The simulations contain proton and iron primaries with energies 18.4 ≤ log 10 (E/eV) ≤ 20.2 in log 10 (E/eV) = 0.2 steps, zenith angles from 65 • to 85 • with a step size of 2.5 • , and 8 equidistantly spaced azimuth angles φ, i.e. coming from geomagnetic East (φ = 0 • ), North-East (φ = 45 • ), North (φ = 90 • ), etc. For each simulation the simulated pulses are located on a star-shaped grid along the v × B and v × (v × B) axes and their bisections, i.e. on concentric rings in the shower plane projected onto the ground. In the following we will refer to the pulses with the same polar angle in the shower plane as "arms" of the star-shaped grid. The spacing is denser close to the shower axis to sample the energy fluence distribution in the 30 to 80 MHz band precisely, and sparser outside, cf. Figs. 1, 2.
The simulation settings match the ambient conditions of the Pierre Auger Observatory [19] which will measure the radio emission from inclined air showers after the upcoming deployment of a large-scale radio detector [4]. The ground plane is set to an observation level of 1400 m above sea level. The atmospheric model fits the average conditions of Malargüe (location of the Pierre Auger Observatory) in October [20]. The refractivity at sea level is set to n 0 − 1 = 3.12 × 10 −4 . We use QGSJetII-04 [21] and UrQMD [22] as high-and lowenergy hadronic interaction models and set a thinning level of 5 × 10 −6 with optimized weight limitation [23].

Fitting the Cherenkov ring
So far there is no established LDF for horizontal air showers that can be used to fit the radio core. New models are currently being developed, e.g. [16], but the results are still being validated. We therefore use a purely geometrical approach that exploits the Cherenkov compression of the radio signal for the core estimation. At a certain angle, the Cherenkov angle, a large fraction of the emission, released during the complete shower evolution, arrives at the same time at ground and enhances the signal strength on a ring around the shower axis, the so-called Cherenkov ring. We here estimate the radio core by a fit of a ring to this feature and define the center of this ring as our radio core. This approach also yields an estimator for the radius of the Cherenkov ring. Its radius depends on the geometrical distance to the source region and on the refractivity in this region [8]. Since geomagnetic and charge-excess emission were found to originate from slightly different regions in the atmosphere [2,24] it is expected that both emission contributions exhibit an independent Cherenkov ring. Hence, in the following we will describe the radio-emission footprint in terms of the geomagnetic energy fluence f geo and charge-excess energy fluence f ce separately. Both can be calculated from the energy fluence in the v × B and v × (v × B) polarisation f v×B , f v×(v×B) for a given core position via (derived from [15]) (1) where Φ denotes the polar angle of the pulse position in the shower plane with respect to the positive v × B axis counting counterclockwise. Equations (1) and (2) are only valid if the pulses of geomagnetic and chargeexcess emission arrive at ground almost simultaneously, i.e., without a significant phase shift giving rise to a circularly polarized signal component. Following the analysis presented in [24] we find an average time delay of 1 ns for pulses within and around the Cherenkov ring. Given pulse widths of tens of nanoseconds this delay is negligible independent of the relative strength of the charge-excess emission or the considered frequency band 1 . Therefore, Eqs. (1) and (2) are valid within the scope of our analysis. In Fig. 2 an example shower with a small geomagnetic angle α (angle between the magnetic field axis and shower axis) is shown. For such showers with a weak geomagnetic emission the interference between geomagnetic and charge-excess emission impacts the position of the maximal fluence and can even completely suppress the Cherenkov ring in the negative v × B half-plane. Note that since the amplitudes of the electric field traces are interfering, the resulting asymmetry in energy fluence, e.g., the squared sum of the amplitudes, is accentuated. For the total fluence no ring can be estimated for signals with Φ = 135 • . In contrast, the geomagnetic and chargeexcess fluences individually exhibit a clear maximum. Note that for Fig. 2 the MC core was used for the calculation of the geomagnetic fluence which does not describe the true core as we will see later. The LDF described by f ce exhibits a broader Cherenkov ring than f geo which motivates to describe both features independently.
Following Eqs. (1) and (2) the calculation of f geo and f ce depends on the core location via the polar angle Φ. Therefore, we fit the ring in an iterative process recalculating f geo and f ce in each step of the minimization. The Cherenkov ring is described by the position along each arm of the star-shaped grid, for which the fluence is maximal. These positions are found using a cubic spline interpolation along each arm of the star-shaped grid. For the minimization process we employ a least squares method with equal weights for each ring position. The calculation of f geo and f ce following Eq. (1) and (2) becomes nonphysical for small values of sin(Φ), hence the v × B axis is excluded. We note that pulses along this axis will not remain at Φ = 0 • respectively Φ = 180 • for the fitted core position and therefore the above equations could provide reasonable energy fluences using the true radio core position. However, a varying number of data points during the fit could result in a bias.
An example fit to the geomagnetic emission for an event with a zenith angle of 85 • coming from North- Fig. 3 Result of the the iterative fit procedure to estimate the radio core. The geomagnetic energy fluence is normalized to the maximum along each arm. For the fit, only the position of the Cherenkov ring along the arm, and not its signal strength, is used. Top: 2D visualization of the fitted Cherenkov ring. For illustration purposes the background constitutes the cubic interpolation of the geomagnetic energy fluence from signals around the Cherenkov ring (signals on or close to the v × B-axis are recovered using the found radio core.). In red the fitted Cherenkov ring and its center, the radio core, are shown. The black star marks the position of the MC core, grey dots show the positions of the simulated pulses. The positions of maximal geomagnetic energy fluence found for each arm of the star-shaped grid are denoted by the black diamonds. Bottom: 1D lateral distribution of the geomagnetic energy fluence for each polar angle of the star-shaped grid except for the v × B axis. Colored points denoted the calculated geomagnetic energy fluence for the simulated pulses. Their interpolation is shown by the dashed lines for each arm, the position of the maximum geomagnetic energy fluence is marked by the colored vertical line. The blue line and box denote the fitted radius of the Cherenkov ring and its uncertainty. The axis distances displayed on the x-axis are calculated using the fitted radio core.
West is shown in Fig. 3. The impact of the underlying interpolation function and the spacing of interpolated points used to find the maxima on each arm is found to be negligible for the obtained results. The displacement between the radio and MC core is estimated to be 125 m in the shower plane. This is a small effect compared to the fitted Cherenkov ring radius of 1198 m, however, due to the high inclination this corresponds to a displacement of 1428 m on ground. The maximal difference between the Cherenkov radii, found on the individual arms, amounts to 40 m. The fit yields an uncertainty of the core displacement in the shower plane of 21 m.
Having two different Cherenkov rings, i.e., the ring in f geo and f ce , encoded in the total signal of f v×B and f v×(v×B) with a similar strength, i.e., for showers with a small sin α, makes it challenging to disentangle them. Hence, in the following we will fit the Cherenkov ring and evaluate a displacement of the emission footprint for the dominating geomagnetic emission only for showers with a larger geomagnetic angle. For a subset of showers with a small geomagnetic angle we will compare the Cherenkov radii independently for both emission contributions.

Investigation of the core displacement for showers with a large sin α
For showers with a large geomagnetic angle, fulfilling sin α > 0.25, we determine the radio core displacement by a fit to the Cherenkov ring. This condition excludes 120 showers coming from North with a zenith angle below 70 • , which we will discuss separately in subsection 3.3. In total 4185 fits yield an accurate result and are analyzed in the following.
We interpret our results as function of the geometric distance d max from the MC core to X max , given by The atmospheric slant depth measured along the shower axis of the ground plane is denoted by X ground , ρ( ) denotes the atmospheric density at the distance along the MC shower axis in the direction of X max . For inclined showers the integral can only be solved numerically as the atmospheric curvature needs to be taken into account. In the first order d max scales with the zenith angle, and only in second order with X max . For a displaced radio core, the geometrical distance between this core and the shower maximum is smaller. However the deviation is of the order of 1 % and therefore negligible for our purposes.
In Fig. 4 we summarize the observed displacement between MC and radio core. In the ground plane we find a displacement of more than 1500 m for the highly inclined showers (Top). This is of the same order as the spacing of the detector stations for the radio upgrade of the Pierre Auger Observatory [4]. To put the magnitude of the core shift into context we also express the offset in the shower plane as a fraction of the fitted radius of the Cherenkov ring which can go up to 15 % (Bottom). The presented core displacement exhibits a pronounced scatter. The sine of the geomagnetic angle, denoted by the color, shows that the found displacement depends on the shower arrival direction. The displacement is strongest for showers coming from West and weakest for East (given the inclination of the magnetic field of ∼36 • both directions translate to the highest sin(α) values). This dependency is further investigated in Fig. 5 where we show the position of the fitted radio core on ground with respect to the MC core in the coordinate origin. We observe a displacement in the direction from which the shower is incoming, i.e., a displacement towards the shower maximum with a small rotation. The previously described scatter manifests as an East-West asymmetry. As the atmosphere in CoREAS simulations is rotationally symmetric these asymmetries in the core displacement cannot be cause by atmospheric properties. We therefore interpret this as a second-order effect introduced by the Earth's magnetic field.

Investigation of the Cherenkov radius in the geomagnetic and charge-excess emission for showers with a small sin α
For simulations with a small geomagnetic angle, i.e. sin α < 0.25, the charge-excess contribution is relatively strong. This allows for an independent analysis of the Cherenkov radius from the charge-excess and geomagnetic emission. As mentioned earlier, for those showers it is not possible to accurately determine the core. Therefore, the core is fixed to the MC core and just the ring radius is fitted with the method given above. Already in Fig. 2 a difference in the Cherenkov radii found for the geomagnetic and charge-excess contributions is visible. This is confirmed for the 120 showers with sin α < 0.25 as shown in Fig. 6. As the core is fixed, all showers are fitted successfully and no further selection is applied. The estimated radius of the Cherenkov ring is systematically larger for the chargeexcess contribution than for the geomagnetic emission. Assuming that the radio emission originates from a point source at an altitude h, the Cherenkov angle is given by θ Ch = arccos(1/n(h)). Hence, one can estimate the height h of the emission region for the given Cherenkov radius. We find that the charge-excess originates from higher up in the atmosphere than the geomagnetic emission. This is in agreement with the results found in [2].
A rotational asymmetry in the charge-excess emission was reported in [2] for one example event. Here, we investigate the charge-excess emission of 120 showers with low sin α and compare the energy fluence found on the Cherenkov ring for all arms of the star-shaped grid (except of the v × B axis) normalized to the average energy fluence over all arms on an event-to-event basis. We find deviations from rotational symmetry with a standard deviation of 9 % in energy fluence. We do not find a convincing proof of a preferred orientation of this asymmetry, but due to the low event number we can also not exclude one. A random spread could possibly be introduced by air-shower sub-structure originating in the hadronic cascade. We note that indications for a possible random deviation from rotational symmetry in the charge-excess emission of inclined air showers have been noticed before 2 .

Comparison of the radio core displacement for different frequency bands
So far we have analysed the radio-emission for frequencies in the 30 to 80 MHz band. This band is used by most current-generation radio experiments and in particular also by the upcoming large-scale Auger radio detector [4]. We now determine the core displacement for footprints in the 50 to 200 MHz frequency band. This is the Fig. 7 Top: Core displacement of the radio-emission footprint at ground for two different frequency bands. Bottom: Fitted radius of the Cherenkov ring in the geomagnetic emission for the two different frequency bands. target frequency band of the GRAND experiment [12], which is also focused on radio measurements of inclined air showers. In Fig. 7 the fitted core displacement (Top) and Cherenkov radius (Bottom) are shown for the two frequency bands.
We find no difference in the average behaviour of the core displacement between the two frequency bands. The spread is smaller for higher frequencies as the Cherenkov ring is more pronounced and thus easier to fit. On average the fitted Cherenkov radius is ∼5 % larger for the higher frequency band. This might be caused by differences in the geometrical distribution of the shower particles which primarily contribute to the radio signal in the considered frequency ranges. However, the exact origin of this deviation needs future investigations which are beyond the scope of this paper.

Interpretation of the displacement as due to refraction
We have shown that for simulations in an atmosphere with a varying refractive index the radio core is systematically displaced from the MC core. In this section we show that this displacement is in agreement with refraction of radio waves in a refractive atmosphere as described by Snell's law (cf. Eq. (4)). For this purpose we develop a model simulating the propagation of a single electromagnetic wave. Furthermore we summarize and validate the treatment of the refractivity in CoREAS and discuss the validity of our model.

Description of refraction using Snell's law
We study the propagation of a single electromagnetic wave through a realistic atmosphere 3 described by a curved trajectory undergoing refraction according to Snell's law. For this purpose we assume discrete changes of the refractive index along the edges of imaginary layers throughout the atmosphere. The propagation within a layer with an upper edge height h i is described by a straight uniform expansion with the phase velocity c n = c 0 /n(h i ) given the refractive index as function of the height above sea level n(h). We adopt the refractive index as frequency-independent (i.e., non-dispersive) in the band from 30-200 MHz that we consider here. The change of direction between two layers with n 1 and n 2 is described in terms of the incidence angle (ϑ) from ϑ 1 to ϑ 2 following Snell's law The refraction is calculated in a curved atmosphere. The relationship between the geometrical distance from ground d g and height above ground h g is given for a zenith angle θ, observation level h obs and the Earth's radius r earth = 6371 km by By solving this quadratic equation, one can calculate the height above sea level for every given distance to the ground by h = h g (d g , θ) + h obs . The refractivity N ≡ n − 1 at a given height is calculated according to the density profile of the given atmospheric model and the given refractivity at sea level (N 0 ), The thickness of each layer is set to 1 m assuring a high accuracy of the calculation 4 .
To predict the magnitude of a core displacement by refraction, we simulate the propagation of an electromagnetic wave along a bent trajectory with an initial direction aligned to the MC axis for a shower with a given zenith angle, towards the ground plane. The intersection between the bent trajectory and the ground plane is compared to the intersection between the ground plane and the MC axis. Given these two points, the core displacement can be inferred as depicted in Fig. 8. In Fig. 9 the predicted core displacement along the ground plane is shown as a function of the geometrical distance along the MC axis for shower geometries with a zenith angle between 65 • and 85 • . The red line symbolises the displacement for a source at a fixed slant depth of 750 g/cm 2 (e.g., shower maximum, average depth of maximum of our set of simulated showers). For a given slant depth, this distance translates to a zenith angle (top x-axis). Our model predicts a core displacement of the order of 1.5 km for the most inclined showers at θ = 85 • . In red squares the displacement is shown for different slant depths between 620 g/cm 2 and 1000 g/cm 2 (typical range in our set of simulated showers) along the MC axis for 5 different zenith angles (θ = 65 • , 75 • , 80 • , 82.5 • , 85 • ). The model predictions are compared to the core displacements determined from the CoREAS simulation set (colored circles: cf. Sec. 3, Fig. 4). The Fig. 9 Comparison between model-predicted and CoREASderived core displacement. Displacement is expressed within the ground plane. The red line constitutes our model prediction for a source at a fixed slant depth of Xmax = 750 g/cm 2 (translation to zenith angle on the top x-axis). The red squares show the displacement as a function of the source slant depth (e.g. Xmax) for each given zenith angle (cf. top x-axis). The colored circles show the displacement determined from the CoREAS simulation set. The color code shows the sine of the geomagnetic angle. displacement is reasonably described by our model in terms of the overall magnitude (red line) as well as the slope as function of the source's slant depth (red squares). Furthermore, our model predicts a displacement always towards the shower incoming direction. This corresponds to a refraction towards the ground, i.e., decreasing angle of incidence, which is given by a radially symmetric atmosphere (cf. Fig. 8). In Fig. 5 this behaviour was also observed for CoREAS simulations as the simulated showers exhibit a radio core displacement almost entirely in the incoming direction of the shower. As emphasised earlier, an East-West asymmetry as seen in CoREAS simulations, cf. Fig. 4 and 5, cannot be described by refractivity.

Refraction and its treatment in CoREAS
For the numerical calculation of the radio emission of an extensive air shower for an observer at ground, the refractive index has to be taken into account for two processes: first, in the generation of the radio emission for each particle; and second, in the the propagation of each electromagnetic wave from a source to an observer. In CoREAS, the former is realistically included in the calculation of the radio emission from each particle using the endpoint formalism [17,27,28]. However, the treatment of the propagation is approximated. Since electromagnetic waves in the radio regime do not suffer from any significant attenuation effects while propagating through air, this propagation is described entirely by two quantities. First, the geometrical distance (d), that the radio wave passes between source and observer, as the intensity of the emission scales with this distance. And second, the light propagation time (t n ) between source and observer which is of crucial importance as it governs the coherence of the signal seen by an observer from the full air shower. In CoREAS, t n is calculated taking into account a refractive index dependent (phase-) velocity of the emission To calculate both quantities, CoREAS assumes a straight path between source and observer (cf. Fig. 8: dashed line). This approximation has implications for the geometrical distance between source and observer as a straight line underestimates the real distance along a curved trajectory. For the calculation of the light propagation time an additional implication arises from the fact that the average refractivity along a straight line varies from the refractivity along a curved trajectory. We find that the average refractivity and consequently the light propagation time is overestimated along straight trajectories.
We stress that the description of the propagation of the radio emission along straight trajectories in CoREAS is not in contradiction with the above-established refraction of the radio emission and the resulting displacement of the radio core in CoREAS simulations. In fact, the refraction of radio waves is a consequence of the fact that the propagation velocity changes with the refractive index c n . It is in fact possible to achieve a displacement of the whole coherent signal pattern at ground by an accurate description of the light propagation time along straight trajectories, as we will demonstrate below.
To verify if the calculation along straight trajectories between source and observer is sufficiently accurate to calculate the radio emission seen from a full extensive air shower, we determine the geometrical distance and light propagation time following bent and straight trajectories (cf. Fig. 8) for several geometries. We simulate the propagation of an electromagnetic wave given incoming direction and atmospheric depth along a bent trajectory towards the ground plane. Once the trajectory intersects with this ground plane the process is stopped and d curved and t curved are calculated via a sum of d i , t n i over all layers. Given the intersection and the initial starting point in the atmosphere, a straight trajectory is defined and d straight and t straight are calculated for comparison.
In Fig. 10 (Top) the geometrical distance is compared between curved and straight trajectories in absolute terms of d curved − d straight , given the ambient conditions used for the above introduced simulation set. The comparison is shown as function of the geometrical distance along the straight trajectory between source and observer. The source positions are set to be at an atmospheric depth of X = 750 g/cm 2 for incoming directions with zenith angles between 65 • and 85 • . We obtain a maximal error of around 4 cm for the most inclined geometries with a path distance of ∼1500 km. With a relative deviation of less than 1 × 10 −6 this approximation is therefore completely suitable.
For the light propagation time, the relative difference for two source positions and one observer position between curved and straight trajectories σ t = ∆t curved − ∆t straight is of relevance as it governs the coherence of the total signal seen by a given observer. We do not have an analytic description for curved trajectories, however we can employ our model to determine the observer position at ground O(P,θ) for every given source position P and initial directionθ. Hence, to find two sources connected with curved trajectories to one observer at ground we have to find the initial direction θ 2 for a given second source which defines a trajectory that connects this source to an observer given by the first source and direction O 1 (P 1 ,θ 1 ). For this purpose we employ a root-finding algorithm that solves the following equation forθ 2 : O 2 (P 2 ,θ 2 ) − O 1 (P 1 ,θ 1 ) = 0 (P 1 , P 2 ,θ 1 are fixed). When the correctθ 2 is found, i.e., O = O 1 = O 2 , the propagation between P 1 or P 2 and O is evaluated for curved and straight trajectories and σ t is calculated. Fig. 10 (Bottom) shows σ t for different geometries and configurations of P 1 and P 2 . The timing error σ t between two sources which are located on one axis with a zenith angle between 65 • and 85 • and depths of 1000 g/cm 2 and 400 g/cm 2 is shown by the blue line. This range of atmospheric depths covers the bulk of the radiation energy release from the longitudinal development of an extensive air shower [2,Fig. 5]. In the most extreme case, the error in the relative arrival times for a source at the beginning of the shower evolution and a source at the end of the shower evolution, estimated using straight tracks, amounts to σ t 0.1 ns. This is well below the oscillation time of electromagnetic waves in the MHz regime. The orange line in the same figure demonstrates the errors made for sources shifted by ±500 m above and below the shower axis along an axis perpendicular to the shower axis at a depth of 750 g/cm 2 . The errors due to the straight-line approximation are even much smaller.
While it may seem paradoxical on first sight that a calculation approximating propagation of electromagnetic waves along straight tracks can yield refractive ray bending, we have shown that the relevant calculation of relative arrival times is described well within the needed Fig. 10 Top: Difference in geometrical distance between a source at a depth of X = 750 g/cm 2 and an observer at ground level for a straight-track calculation and curved-track calculation for zenith angles from 65 • to 85 • . Bottom: Difference in the relative arrival times calculated with curved and straight-track propagation σt = ∆t curved − ∆t straight arising for two source positions and one observer position. Distance to Xmax is calculated using a fixed depth of Xmax = 750 g/cm 2 . accuracy, i.e., is fully adequate for this purpose. We note that, similar to our findings, it was already found based on analytic calculations in reference [29, cf. Fig  9] that a straight-line approximation is sufficient for the calculation of relative arrival times of radio waves in extensive air showers.
Additionally the refractive ray bending changes the incoming direction of the radio emission. This has implications for the reconstruction of the radio emission with real radio antennas as their response pattern is direction-dependent. We find a maximum change of direction of ∼0.14 • , which is in agreement with [30]. This is below current experimental accuracy as well as the change in the incoming direction between early and late observers on the ground plane, estimated as O(1 • ) for a 85 • shower.

Conclusions
We have established that the radio core is displaced with respect to the MC core in CoREAS simulations of inclined air showers. This displacement shows no significant dependence on the considered frequency band for non-dispersive refractivity. We have developed a model which reproduces this displacement quantitatively, describing it as a result of refraction of the radio waves during propagation in an atmosphere with a refractive index gradient. We have also discussed the validity of approximations made in CoREAS and shown that they are adequate to describe refractive effects in air-shower radio simulations. We have found indications that there are secondary effects causing core displacements which are not related to the refractive index of the atmosphere or the atmosphere in principle but rather must be related to the geomagnetic field.
These findings have several implications towards the observation of cosmic rays with the radio detection technique. For the development of reconstruction algorithms, assuming the MC core as symmetry center of the radioemission footprint will disturb its lateral distribution, causing a mismodelling of the signal distribution. In observations the refractive displacement primarily has to be taken into account in the interpretation of the reconstructed shower geometry. Considering hybrid detection and reconstruction, refractive displacement has to be taken into account when comparing/combining results across different detection techniques.