The effect of the atmospheric refractive index on the radio signal of extensive air showers

TL;DR

This study quantifies how natural variations in atmospheric refractive index affect radio measurements of air shower properties, revealing systematic errors in the inferred depth of shower maximum that can be mitigated for improved accuracy.

Contribution

The paper evaluates the impact of refractive index variations on radio measurements of air showers using simulations, providing correction estimates for more precise $X_{max}$ determination.

Findings

A 4% change in refractivity causes a 3.5-11 g/cm^2 error in $X_{max}$ at 30-80 MHz.

At 120-250 MHz, the error increases to 10-22 g/cm^2.

A simple Cherenkov angle model qualitatively agrees with simulation results.

Abstract

For the interpretation of measurements of radio emission from extensive air showers, an important systematic uncertainty arises from natural variations of the atmospheric refractive index $n$. At a given altitude, the refractivity $N=10^6\, (n-1)$ can have relative variations on the order of $10 \%$ depending on temperature, humidity, and air pressure. Typical corrections to be applied to $N$ are about $4\%$. Using CoREAS simulations of radio emission from air showers, we have evaluated the effect of varying $N$ on measurements of the depth of shower maximum $X_{\rm max}$. For an observation band of 30 to 80 MHz, a difference of $4 \%$ in refractivity gives rise to a systematic error in the inferred $X_{\rm max}$ between 3.5 and 11 $\mathrm{g/cm^2}$, for proton showers with zenith angles ranging from 15 to 50 degrees. At higher frequencies, from 120 to 250 MHz, the offset ranges from 10 to 22 $\mathrm{g/cm^2}$. These offsets were found to be proportional to the geometric distance to $X_{\rm max}$. We have compared the results to a simple model based on the Cherenkov angle. For the 120 to 250 MHz band, the model is in qualitative agreement with the simulations. In typical circumstances, we find a slight decrease in $X_{\rm max}$ compared to the default refractivity treatment in CoREAS. While this is within commonly treated systematic uncertainties, accounting for it explicitly improves the accuracy of $X_{\rm max}$ measurements.

Equations

$$\cos \alpha$$

$$\alpha$$

$$X_{\max }\equiv X(h_0)=\frac{1}{\cos \theta }{\int }_{h_0}^{\infty }\rho (h)\mathrm{d}h,$$

$$X(h_0)=\frac{10}{g}\frac{p(h_0)}{\cos \theta },$$

$$\sqrt{N_1/N_0}{h}_1\alpha (h_1)=h_0\alpha (h_0),$$

$$T(h)$$

$$p(h)$$

$$p_{\mathrm{sat}}$$

$$p_w$$

$$N=77.6890\frac{p_d}{T}+71.2952\frac{p_w}{T}+375463\frac{p_w}{T^2},$$

$$\mathrm{MSE}(k)=\frac{1}{N_{\mathrm{ant}}}\sum _{\mathrm{antennas}}{\left(A{f}_{\mathrm{test}}(x,y)-f_k(x,y)\right)}^2,$$

$$A=\frac{\sum f_{\mathrm{test}}(x,y)f_k(x,y)}{\sum f_{\mathrm{test}}(x,y{)}^2}.$$

$$\Delta X_{10\%}=3.00\left(\frac{R}{1\mathrm{km}}\right)-1.37g/c{\mathrm{m}}^2,$$

$$\Delta X_{10\%}=4.18\left(\frac{R}{1\mathrm{km}}\right)+11.24g/c{\mathrm{m}}^2.$$

Full text

††journal: Astroparticle Physics

1 Introduction

In recent years, the techniques for measuring and modelling radio emission from air showers induced by high-energy cosmic rays have developed rapidly [1]. The radio detection method has achieved high precision in estimating the air shower and primary particle properties [2, 3, 4] which allows for very precise measurements of the primary particle type and the energy of the primary cosmic ray [5]. In particular, the LOFAR radio telescope [6] has been used successfully for this, due to its densely instrumented core region located in the Netherlands. In an inner ring of $320\,\mathrm{m}$ diameter, we can use 288 low-band antennas, measuring in the $30-80\,\mathrm{MHz}$ range, for cosmic-ray measurements. Additionally, there are also 288 high-band antennas measuring in the $110$ to $190\,\mathrm{MHz}$ range. In an extended core region of about $6\,\mathrm{km^{2}}$, nearly 1800 additional low-band antennas have been installed, grouped into stations of 96. Depending on strength and location of the air shower signal, up to four stations outside the inner ring can be used as well. The signals from air showers have been routinely measured with LOFAR since 2011 [7].

In an air shower, secondary electrons and positrons are produced, which undergo charge separation as they travel through the Earth’s magnetic field. This leads to transverse currents producing radio emission. This emission reaches the ground as a short pulse on the order of 10 to $100\,\mathrm{ns}$ long, with a specific lateral intensity distribution or ‘footprint’ [8] that depends on the depth of shower maximum $X_{\rm max}$. The number of produced secondary particles peaks at $X_{\rm max}$. This point, expressed as the column density of traversed matter ($\mathrm{g/cm^{2}}$), varies with primary particle type and is therefore an important quantity to measure in composition studies.

The measured lateral intensity distributions are compared to microscopic Monte Carlo simulations of air showers, to infer the properties of the primary cosmic ray. To simulate the radio signal at the antennas we use CoREAS [9], a simulation of the radio emission from the individual particles in the cascade simulated with CORSIKA [10]. Fitting these simulated radio footprints to measured air showers allows us to infer $X_{\rm max}$ to a precision of $20\,\mathrm{g/cm^{2}}$. As this precision is comparable to that of fluorescence detection [11, 12], it is well suited for composition studies. The radio detection method is therefore a complementary technique, as it is not limited to dark clear nights, and because its duty cycle is limited only by technical conditions and thunderstorms, it can reach in principle up to $100\,\mathrm{\%}$.

The detected radio signal depends on the difference in travel times of radio waves and particles. Therefore, it is important to apply an accurate value of the refractive index $n$. Variations in the refractive index lead to changes in the radio intensity footprint on the ground because the angle of peak emission depends on $n$. The refractive index, which at sea level is about $n\approx 1.0003$, exhibits natural variations due to weather conditions, at the level of the fifth decimal. It is therefore common to define the refractivity $N=10^{6}\,(n-1)$, which emphasizes relative variations as these depend on $(n-1)$.

In this paper, we quantify the influence of the refractive index variations on the depth of shower maximum, to reduce the systematic uncertainty of the $X_{\rm max}$ measurements. In particular, CoREAS currently assumes a constant default value of the refractive index at each altitude in the atmosphere, and we explore how much the systematic error can be reduced by a more accurate treatment of the local atmosphere. In the next section, we describe a toy model for the radiation from air showers, which qualitatively explains how unaccounted variations in the refractive index give rise to uncertainties in determining $X_{\rm max}$. In Sect. 3 we review the equations used to describe atmospheric parameters and their altitude profiles, and to calculate the refractive index. In Sect. 4, the method of fitting intensity distributions is described, and Sect. 5 gives the results for the systematic offsets of $X_{\rm max}$.

2 Toy model for the effect of the refractive index on radiation from air showers

The depth of shower maximum $X_{\rm max}$ can be inferred from the radio intensity footprint measured on the ground. In this section we show the radio footprint changes with the Cherenkov angle, which is a function of the refractive index. In particular, if the refractive index is higher than expected, the method based on the intensity footprint will underestimate $X_{\rm max}$.

In an extensive air shower, the magnetic field of the Earth induces an electric current, as the Lorentz force has opposite direction for electrons and positrons in the shower front. This current is transverse to the direction of the shower. The number of electrons and positrons depends strongly on the interaction depth of the shower and peaks at $X_{\rm max}$. The induced current is therefore strongly time dependent and emits electromagnetic radiation. The shower front has a thickness on the order of meters, and the radiation is coherent at wavelengths longer than this thickness, i.e. at radio wavelengths [13]. The transverse current resides in the shower front and thus moves towards the Earth surface with a velocity exceeding the speed of light in air [14]. Therefore, radio waves are emitted because of coherent Cherenkov emission [15]. In addition to the emission from the transverse current there is also a smaller contribution from the net negative charge buildup in the shower front. From the polarization of the radio signal, this contribution was found to be $11\,\mathrm{\%}$ on average at the LOFAR site [16].

In the $30-80\,\mathrm{MHz}$ band primarily used at LOFAR, the emission along the Cherenkov angle and the non-Cherenkov emission have roughly the same magnitude, while at higher frequencies, the Cherenkov emission dominates (see Fig. 4 in [17]). This has been confirmed by LOFAR observations above $110\,\mathrm{MHz}$ where a clear ring-like emission pattern is found [18]. The radius of the ring was on the order of $100\,\mathrm{m}$ as is expected from the Cherenkov angle and the distance to $X_{\rm max}$. It has also been observed at GHz frequencies by the CROME [19] and ANITA [20] experiments.

The angle $\alpha$ where the Cherenkov emission peaks, is given by

[TABLE]

where $\beta=v/c$ is the velocity of the shower front with $\beta=1$ to a good approximation, the refractive index of air $n\approx 1.0003$ at sea level, and it varies with altitude. For convenience we also use the refractivity $N$ throughout the text, as we will consider relative variations in $N$.

The depth of shower maximum $X_{\rm max}$ is reached at an altitude $h_{0}$, which is given by the altitude-dependent density $\rho(h)$ and the zenith angle $\theta$. The relation is [21]

[TABLE]

where $X$ is the column density expressed in $\mathrm{g/cm^{2}}$. Therefore, $X$ is also proportional to the pressure,

[TABLE]

with $g$ the gravitational acceleration and $p$ the pressure in $\mathrm{Pa}$.

As a toy model for analyzing the effect of varying refractive index, we use the approximation that the size of the radio footprint is proportional to the base of a cone located in the shower plane, with a half-opening angle of $\alpha$ with respect to the the direction of the incoming primary particle. The shower plane is defined as the plane perpendicular to the incoming direction of the primary particle. Moreover, we assume that all radiation is produced near $X_{\rm max}$. As a consequence, the size of the radio footprint on the ground would be proportional to the geometric distance to $X_{\rm max}$ and to the Cherenkov angle at the altitude of $X_{\rm max}$. Variations in refractive index $n$ at altitude $h_{0}$ would then translate to variations in radio footprint size via the Cherenkov angle. The radio footprint with its non-circular symmetric structure [22, 8] falls off smoothly with distance, hence it has no sharply defined “size”. However, the important point for this model is that a given footprint would scale, both with distance to $X_{\rm max}$ and with the Cherenkov angle.

The model is expected to be more accurate at frequencies above $100\,\mathrm{MHz}$ where the Cherenkov mechanism dominates. At lower frequencies, the intensity pattern depends less strongly on the Cherenkov angle, and therefore on the refractive index.

As the refractive index is usually discussed in context of optical refraction, it should be noted that the additional effect of variations in $n$ on bending of signal propagation paths (through Snell’s law) is negligible for us. A signal path traveling from a medium with $n=1$ to a medium with $n=1.0003$ at a 60 degree incidence angle will be bent by 0.03 degrees. This is already below the resolution of about 0.1 degree attainable with LOFAR [23], and natural variations in $n$ are still an order of magnitude smaller.

In composition studies, we fit radio intensity profiles simulated with CoREAS to the measured intensity profile on the ground [2], to estimate $X_{\rm max}$. In Fig. 1, the effect of an increase in $n$ is shown schematically, for a proton primary particle of $10^{17}\,\mathrm{eV}$ from zenith. These have an average $X_{\rm max}$ $\approx 670\,\mathrm{g/cm^{2}}$ (from CoREAS simulations), corresponding to an altitude of $3.51\,\mathrm{km}$. If the refractive index is higher than expected, fitting the intensity profiles at ground will underestimate $X_{\rm max}$ (blue lines), as for the actual refractive index, the Cherenkov angle is larger (red lines). The lower altitude in the atmosphere corresponds to a higher $X_{\rm max}$; the modeled difference amounts to $17\,\mathrm{g/cm^{2}}$ for a $10\,\mathrm{\%}$ increase in refractivity $N$. This is therefore a systematic uncertainty on $X_{\rm max}$.

The offset in $X_{\rm max}$ follows from the condition that the footprint size, which in this approximation is just the diameter of the intersection of the Cherenkov cone with the ground or shower plane, is kept constant. This gives

[TABLE]

with $h_{0}$ and $h_{1}$ the altitude of the $X_{\rm max}$ point at fixed footprint size, for standard refractivity $N_{0}$ and $N_{1}=1.1\,N_{0}$ respectively. The square root arises from the small-angle approximation in Eq. 2; we have also taken $\sin\alpha\approx\alpha$. It holds both for vertical and for inclined showers and takes into account that $N$, and therefore $\alpha$, varies with altitude as shown in the equations in the next section.

In summary, the steps for calculating the modelled shift in $X_{\rm max}$ given an original $X_{\rm max}$ at altitude $h_{0}$, are:

establish $N_{0}$ at altitude $h_{0}$ using the equations in the next section (Eqs. 6 through 10).

• 2.

consider a different refractivity profile with altitude, e.g. $N_{1}(h)=1.1\,N_{0}(h)$

• 3.

solve $h_{1}$ from Eq. 5

• 4.

obtain the column density (shifted $X_{\rm max}$) at this height, from Eq. 4 and Eq. 7 below

In the example of Fig. 1, we would have $h_{0}=3.51\,\mathrm{km}$, for which $N_{0}=192.3$. When considering a $10\,\mathrm{\%}$ higher refractivity, i.e. $N_{1}(h)=1.1\,N_{0}(h)$ at all altitudes, we obtain $h_{1}=3.31\,\mathrm{km}$ when keeping the footprint size constant. At this level, we would have a column density $X=687\,\mathrm{g/cm^{2}}$, which is a systematic offset of $17\,\mathrm{g/cm^{2}}$.

3 The atmospheric refractive index

The refractive index $n$ and refractivity $N$ vary with temperature, pressure, and humidity in the atmosphere. As these depend on altitude, we use a parametrization for pressure and temperature corresponding to the US Standard atmosphere [24], where we have rewritten the equations to a slightly more compact form. We use only the bottom layer of this model, from 0 to $h=11\,\mathrm{km}$, where $h=0$ defines sea level. It is valid in the troposphere under ‘average’ circumstances,

[TABLE]

Here, $p_{0}=1013.25\,\mathrm{hPa}$ and $T_{0}=288.15\,\mathrm{K}$ are standard sea level pressure and temperature. These can be varied according to local circumstances. The temperature lapse rate $L$ is assumed constant at $L=6.5\,\mathrm{K/km}$. The remaining constants are the gravitational acceleration $g=9.80665\,\mathrm{m/s^{2}}$, the ideal gas constant $R=8.31447\,\mathrm{J/(mol\,K)}$, and the molar mass of dry air $M=0.0289644\,\mathrm{kg/mol}$. In the limit $L\rightarrow 0$, Eq. 7 reduces to the familiar exponential barometric formula.

Given the relative humidity $H$, the partial pressure of water vapor is calculated using the Magnus formula for the saturation pressure [25]:

[TABLE]

where $t$ is temperature in $\mathrm{{}^{\circ}C}$ (in contrast to $T$ in K), and $p_{\rm sat}$ is the saturation pressure for water vapor. For the empirically determined constants we take the values from [25], which are $a=6.1121\,\mathrm{hPa}$, $b=17.502$, $c=240.97\,\mathrm{K}$ for temperatures above $0\,\mathrm{{}^{\circ}C}$, and $a=6.1115\,\mathrm{hPa}$, $b=22.452$, $c=272.55\,\mathrm{K}$ below $0\,\mathrm{{}^{\circ}C}$. The relative uncertainty in $p_{\rm sat}$ is then given as $0.1\,\mathrm{\%}$ over the range of $-50$ to $+40\,\mathrm{{}^{\circ}C}$ .

The radio refractivity is parametrized according to [26] as

[TABLE]

where $T$ is temperature in K; $p_{d}$ and $p_{w}$ are the partial pressures (hPa) of dry air and of water vapor, respectively. The total air pressure is $p=p_{w}+p_{d}$. The small influence of carbon dioxide is included in the dry-air contribution. The accuracy of this equation is given as $0.02\,\mathrm{\%}$ for the first term in Eq. 10, and $0.2\,\mathrm{\%}$ for the second and third term combined. This evaluates to a relative uncertainty of less than $0.1\,\mathrm{\%}$ in total, which is sufficient for our purposes.

It should be noted that at the radio frequencies of interest here, the refractivity values are different from those at infrared, visible, and UV wavelengths, such as considered in [21] for the fluorescence detection technique at Pierre Auger Observatory. In particular, the presence of water vapor significantly raises the radio refractivity, while it tends to lower the infrared refractivity slightly below that of dry air. The latter is depicted in [27] regarding optical Cherenkov telescopes. Therefore, the precision formulas of Edlén (updated in [28]) and Ciddor [29], defined for visible and near-infrared wavelengths are not applicable for us, and one can also not simply take the infinite-wavelength limit of those. For instance, for air at standard pressure, $20\,\mathrm{{}^{\circ}C}$, and $50\,\mathrm{\%}$ relative humidity, the Ciddor equation gives $N=268$, while Eq. 10 gives $N=319$. It follows that for radio detection, accounting for humidity is more important than for fluorescence detection.

With the above formulas and definitions, we calculate the atmospheric profile of the refractivity versus altitude for different values of temperature, pressure, and humidity. In the US Standard atmosphere one uses the geopotential altitude, which takes into account the decrease in gravitational acceleration $g$ with altitude. As the difference between geometric and geopotential altitudes at $h=5\,\mathrm{km}$ is only 4 m, we do not correct for the difference.

An example plot of the altitude dependence of $N$ is shown in Fig. 2 for two values of temperature and humidity. It is clear that humidity cannot be neglected, especially in the Netherlands where relative humidity near sea level is on average roughly $80\,\mathrm{\%}$ [30].

In CoREAS, the refractivity is set to a default value of $N=292$ at sea level, and scaled with density according to the US Standard atmosphere, which considers dry air and $T=15\,\mathrm{{}^{\circ}C}$ at sea level. In Fig. 3(a), we show the relative correction factor to $N$ with respect to its default value, at an altitude of $3.5\,\mathrm{km}$. This corresponds to $X_{\rm max}$ = $670\,\mathrm{g/cm^{2}}$, which is the average value found for protons of $10^{17}\,\mathrm{eV}$ energy. The correction is plotted against the ground temperature for several values of the relative humidity at the given altitude. At a sea-level temperature of $10\,\mathrm{{}^{\circ}C}$ and $70\,\mathrm{\%}$ relative humidity, the true refractivity is about $2\,\mathrm{\%}$ lower than the default value. At higher sea-level temperatures, humidity plays a larger role, and the curves show a larger spread.

For air showers coming in at a 45 degree zenith angle, the altitude of $X_{\rm max}$ is about $6.1\,\mathrm{km}$. At this altitude, the actual $N$ is around $5\,\mathrm{\%}$ lower than the default value.

For the air pressure, we have taken the standard value of $p=1013.25\,\mathrm{hPa}$. Natural variations, which are on the order of $\pm 2\,\mathrm{\%}$, have the effect of lowering or raising the altitude of a given column density $X_{\rm max}$.

The relative humidity is expected to vary with altitude, dependent on conditions like cloud cover. The humidity in a region around the $X_{\rm max}$ altitude will be the most important. Therefore, to accurately estimate $N$, it is necessary to use atmospheric profile data such as available e.g. through the Global Data Assimilation System (GDAS) [31]. Uncertainties in these data translate into the uncertainty in $N$. Ref. [21] gives a comparison between GDAS data and weather balloon measurements in Argentina. Over the years 2009 and 2010, and in the altitude range of 3.5 to $6\,\mathrm{km}$, this difference is $\pm 0.5\,\mathrm{{}^{\circ}C}$ for temperature, $0.5\,\mathrm{hPa}$ for pressure, and $0.05\,\mathrm{hPa}$ for water vapor pressure. The latter corresponds to about $2\,\mathrm{\%}$ relative humidity at $3.5\,\mathrm{km}$ altitude, and about $7\,\mathrm{\%}$ at $6\,\mathrm{km}$.

This shows that the GDAS data accurately represent the local circumstances. Although the uncertainty numbers may be different for other experiment sites, they are expected to be on the same order. The resulting relative uncertainty in $N$ is around $0.5\,\mathrm{\%}$ and is dominated by the humidity uncertainty. This is sufficient for our purposes.

4 Method

To investigate the effect of changing the refractive index on $X_{\rm max}$ measurements, we considered proton showers with a primary energy of $10^{17}\,\mathrm{eV}$ for five different zenith angles. For each case, we have generated an ensemble of 50 simulated showers with $N$ at its default value of $N=292$ at sea level. We have also generated another set of 50 showers where the refractivity $N$ is higher by $\mathrm{10}{\%}$ at all altitudes. This is done using the same random number seeds, ensuring the evolution of the particles is exactly the same. Only the radio emission is recalculated, taking into account the higher refractivity. This method is similar to the one used in [2] that was used for the composition analysis at LOFAR [5]. The difference is that here we compare two simulated ensembles instead of comparing simulations to measured data.

For all showers, we evaluate the lateral distribution of signal intensity in the shower plane $f(x,y)$, for a star-shaped pattern of antennas. In this pattern, 160 antennas are laid out along 4 lines, at a distance $25\,\mathrm{m}$ apart. The lines intersect at the origin and make angles of 45 degrees to each other, forming an octagonal star pattern. For the signal intensity we use a bandpass filter to limit the frequency range to $30-80\,\mathrm{MHz}$, relevant for the LOFAR cosmic-ray project. For comparison, we also consider the $120-250\,\mathrm{MHz}$ band, where Cherenkov effects are more important.

We then take one of the showers with higher refractivity as a ‘test shower’, and fit the lateral distribution of each of the 49 other showers at the default refractivity to it. When fitting the shower with index $k$ to the test shower, this yields a mean-square error as a fit quality measure:

[TABLE]

where ($x$, $y$) are the antenna positions in the shower plane, and $A$ is a scale factor that is taken as a free parameter. This is proportional to a reduced $\chi^{2}$ for the case where all antennas have the same uncertainty on the intensity. Additionally, as we simulate antennas in a star-shaped pattern, we apply weight factors such that each antenna represents the same amount of area in the footprint. The optimal value of $A$ that minimizes the MSE, is given by

[TABLE]

For every shower, we plot the fit quality (MSE) versus $X_{\rm max}$. It is expected to have a minimum, and to lowest order, to have a quadratic dependence around the minimum. Therefore, we fit the points with a parabola. The location of the minimum of the parabola is used as an estimator for $X_{\rm max}$. For one shower this is shown in Fig. 4, for a limited range around the minimum. The scatter around the fitted parabola arises from natural shower-to-shower fluctuations.

We found it useful to weight the fit-quality data points like shown in Fig. 4 by their inverse square in the parabolic fit, thus putting more emphasis on well-fitting profiles. This lowers the uncertainty in estimating the minimum $X_{\rm max}$ by up to $20\,\mathrm{\%}$, without introducing a bias. Proceeding this way, the inferred $X_{\rm max}$ of the test shower with higher refractive index has a systematic offset with respect to its true $X_{\rm max}$ value. By taking in turn each of the 50 showers in the ensemble as test shower, this offset is evaluated along with its statistical uncertainty. This quantifies the effect of variations in the refractive index in e.g. a composition analysis where simulations are fitted to data.

We have limited the range of the fit to include those values of $X_{\rm max}$ within $\pm 40\,\mathrm{g/cm^{2}}$ of the expected value. As this range is determined by the true $X_{\rm max}$ and the offset, we shift the fit range in a second iteration. The offsets found in both iterations are consistent. We have discarded 5 showers at the low, and 5 at the high end of the true $X_{\rm max}$ range, as these give less accurate parabola fits due to lack of data points at one side of the range.

The parabolic fit procedure has also been tested using only showers with the same refractive index. On average over 50 showers, the true $X_{\rm max}$ is reproduced within standard errors. For the ensemble of 50 minus 10 showers, the standard error ranges from 1 to $2\,\mathrm{g/cm^{2}}$ for zenith angles from 15 to 50 degrees, respectively.

5 Results

We have evaluated the bias in $X_{\rm max}$ for a $10\,\mathrm{\%}$ increase in refractivity, for five different zenith angles between 15 and 50 degrees. An increase in $N$, when not taken into account, leads to an overestimate in the geometric distance to $X_{\rm max}$, and hence to an underestimate in $X_{\rm max}$. Making the simplistic assumption that all radiation is produced at $X_{\rm max}$, and that the size of the radio footprint scales with the Cherenkov angle (Eq. 1), we apply the procedure given in Sect. 2 to find the shift in $X_{\rm max}$.

For the $30-80\,\mathrm{MHz}$ band, the results from the simulation fits are plotted together with the model predictions in Fig. 5(a). A linear fit to the results for individual showers is also shown (blue line). The modeled offset in $X_{\rm max}$ is nearly linear in the distance to $X_{\rm max}$, and is a bit less than twice the offsets from simulations. The modeled offset reproduces the general dependence on geometric distance, but the simulated offsets are about a factor 2 smaller. As explained in Sect. 2, this is understood, as for low frequencies the non-Cherenkov emission is about equally important as the Cherenkov emission, leading to a weaker dependency on $N$. The offset ranges from about $9\,\mathrm{g/cm^{2}}$ for a zenith angle of 15 degrees, to $28\,\mathrm{g/cm^{2}}$ at $\theta=50$ degrees. A linear fit suffices to describe the data; a higher-order polynomial curve does not significantly reduce the residuals. The standard deviation of the residuals is $8.9\,\mathrm{g/cm^{2}}$, which is also the intrinsic uncertainty in this fit method, due to natural shower-to-shower fluctuations. The relation for the shift $\Delta X$, defined as the underestimation of $X_{\rm max}$ per $10\,\mathrm{\%}$ increase in refractivity, is

[TABLE]

with $R$ the geometric distance to $X_{\rm max}$. It can be scaled with the relative change in refractivity.

For the $120-250\,\mathrm{MHz}$ band, the $X_{\rm max}$ offset is shown in Fig. 5(b). Here, the simulated offsets are well reproduced by the model, up to an almost constant additional shift. The fitted linear function for the offset is

[TABLE]

The standard deviation of the fit residuals is $7.2\,\mathrm{g/cm^{2}}$ which is again not significantly reduced when using a quadratic function.

The relative correction to the refractivity from Fig. 3 is almost always below unity for vertical showers, and always below unity for inclined showers in the temperature range relevant at LOFAR. Therefore, this gives a systematic offset, not only in individual $X_{\rm max}$ estimates, but also in the average $X_{\rm max}$ inferred from many air showers.

Over the range $-5$ to $20\,\mathrm{{}^{\circ}C}$, the actual $N$ is about $1.5\,\mathrm{\%}$ below the default simulated value for vertical showers. For inclined showers, the actual $N$ is about $4.5\,\mathrm{\%}$ below the default. Hence, the average $X_{\rm max}$ is overestimated by about $1.5\,\mathrm{g/cm^{2}}$ for near-vertical showers, and by about $11\,\mathrm{g/cm^{2}}$ at $45\,\mathrm{degrees}$ inclination.

The bias is larger for more inclined air showers; it has been noted in the LOFAR composition study [5] that more inclined air showers (above 32 degrees zenith angle) had on average a higher inferred $X_{\rm max}$ than the more vertical ones. The difference in average $X_{\rm max}$ amounts to $16\,\mathrm{g/cm^{2}}$ between the two cases. From the offsets in Eq. 13, it follows that at most $9\,\mathrm{g/cm^{2}}$ of this can be explained by a bias caused by using an incorrect value of the refractive index.

The fact that $N$ is slightly overestimated in simulations is also relevant for other radio detection experiments such as AERA [3] and Tunka-Rex [4]. These experiments, located in Argentina and Siberia respectively, have a temperature range different from the LOFAR site. Nevertheless, for the case of inclined showers, the relative correction as from Fig. 3(b) is below unity down to well below $-40\,\mathrm{{}^{\circ}C}$.

The numerical constants in Eqs. 13 and 14 follow from the simulations, which depend on location parameters such as height above sea level and the geomagnetic field vector. Therefore, for other experiments the bias on the average $X_{\rm max}$ is expected to be on the same order, but with slight differences due to variations in location-specific parameters, leading to different constants in Eqs. 13 and 14.

6 Summary

The technique of measuring radio signals from air showers to infer the mass composition of cosmic rays, relies on accurate measurements of the depth of shower maximum $X_{\rm max}$. One of the systematic uncertainties on $X_{\rm max}$ is given by the refractive index of air in the atmosphere, which exhibits natural variations.

We have evaluated the effect of variations in the refractive index of air on determining $X_{\rm max}$. Using a procedure similar to that used in the composition study at LOFAR [5], we have taken simulated proton showers at a $10\,\mathrm{\%}$ higher refractivity $N=10^{6}\,(n-1)$. We have fitted them with a Monte Carlo ensemble of 49 showers simulated at a default value of $N$, excluding the one corresponding to the higher-$N$ shower being fitted. The minimum in the least-squares fit quality yields, on average over many showers, the systematic offset in $X_{\rm max}$.

These offsets were found to be proportional to the geometric distance to $X_{\rm max}$, for zenith angles ranging from 15 to 50 degrees. The effect is roughly twice as strong for the $120-250\,\mathrm{MHz}$ band as for the $30-80\,\mathrm{MHz}$ band. Given variations in $N$ on the order of $4\,\mathrm{\%}$, from the atmospheric effects described in Sect. 3, the offsets would range from about $3.5$ to $11\,\mathrm{g/cm^{2}}$ for $30-80\,\mathrm{MHz}$, and from $10$ to $22\,\mathrm{g/cm^{2}}$ for $120-250\,\mathrm{MHz}$.

A simplistic model in which the radio intensity pattern is assumed to be proportional in size to the Cherenkov cone starting from the $X_{\rm max}$ point, qualitatively describes the effect; the fitted offsets were found to be just above half the modeled offsets in the $30-80\,\mathrm{MHz}$ band, and about 20 % above the modeled offsets in the $120-250\,\mathrm{MHz}$ band.

Calculated profiles of refractivity versus altitude show that one cannot use a single default value of $N$ in simulations. This leads to a bias in the average $X_{\rm max}$, depending on the choice of a constant $N$ either for near-vertical or for inclined showers. The accuracy can be improved for individual showers by using Eqs. 13 and 14.

As a further improvement, it would be required to include detailed atmospheric data, including the particular refractivity profile at the time of the air shower, into a next version of Corsika / CoREAS. The GDAS database [31] is useful for this. This allows to fully account for the refractive index variations, and to re-evaluate the LOFAR measurements at the best level of detail. This is the subject of a future publication.

Acknowledgements

We acknowledge financial support from the Netherlands Organization for Scientific Research (NWO), VENI grant 639-041-130, the Netherlands Research School for Astronomy (NOVA), the Samenwerkingsverband Noord-Nederland (SNN) and the Foundation for Fundamental Research on Matter (FOM). We acknowledge funding from the European Research Council under the European Union’s Seventh Framework Program (FP/2007-2013) / ERC (grant agreement n. 227610) and under the European Union’s Horizon 2020 research and innovation programme (grant agreement n. 640130).