Uncertainties in the Estimation of Air Shower Observables from Monte Carlo Simulation of Radio Emission

TL;DR

The paper systematically compares two leading Monte Carlo codes, CoREAS and ZHAireS, for radio emission from extensive air showers by using identical primaries, geometries, magnetic fields, and a consistent Gladstone-Dale refractive-index model. It evaluates electric-field observables, energy radiated in the radio band, and the reconstruction of the depth of shower maximum, employing a strategy to minimize shower-to-shower fluctuations. The results show strong agreement between the codes across most observables, with typical electric-field differences around a few percent, radiation-energy differences at the few-percent level, and $X_{ m max}$ reconstruction differences of about $~{ m 10}$ g cm$^{-2}$. Differences at higher frequencies and near the Cherenkov ring become larger due to coherence loss and code-specific thinning; overall, the study validates the reliability of radio-based EAS analyses while highlighting areas for further investigation. The authors advocate cross-validation with both CoREAS and ZHAireS in data analyses to ensure robust inference of primary energy and mass composition.

Abstract

The detection of extensive air showers (EAS) induced by cosmic rays via radio signals has undergone significant advancements in the last two decades. Numerous ultra-high energy cosmic ray experiments routinely capture radio pulses in the MHz to GHz frequency range emitted by EAS. The Monte Carlo simulation of these radio pulses is crucial to enable an accurate reconstruction of the primary cosmic ray energy and to infer the composition of the primary particles. In this work, a comprehensive comparison of the predicted electric field in EAS simulated with CoREAS and ZHAireS was conducted to estimate the systematic uncertainties arising from the use of different simulation packages in the determination of two key shower observables namely, the electromagnetic energy of the EAS and the depth of maximum development ($X_{\rm max}$). For this comparison, input parameters and settings as similar as possible were used in both simulations, along with the same realistic atmospheric refractive index depending on altitude, which is crucial for the prediction of radio emission properties of EAS. In addition, simulated EAS with very similar values of depth of maximum development were selected. Good agreement was found between CoREAS and ZHAireS, with discrepancies in the dominant electric field components generally remaining below 10\% across the frequency range of a few MHz to hundreds of MHz, relevant for most radio detection experiments, translating into uncertainties in the determination of energy below $5\%$ and $\simeq 10\,\mathrm{g/cm^2}$ in $X_{\rm max}$. Our work underscores the need for further studies to clarify their origin and impact on $X_{\rm max}$ inference in composition analyses.

Figures (10)

• Figure 1: Sketch of the coordinate system used in this work defining north, south, east, and west in relation to the Earth's magnetic field orientation. The configurations of the magnetic field (horizontal or vertical) adopted in this study are displayed (green and light blue arrows). The antenna positions typically along the north-south ($x$--axis) and east-west ($y$--axis) direction, unless otherwise indicated, are also shown. The shower axis, by default along the north-south direction for showers arriving at the ground from north, and the zenith angle $\theta$ are indicated. The curvature of the Earth, not shown in the sketch, is accounted for in the simulations. The Cherenkov ring of the radio signal, defined in the text, is represented as a blue ellipse.
• Figure 3: Modulus of the Fourier component $E_y$ of the electric field as obtained in simulations performed with ZHAireS and CoREAS, 20 showers each, for primary protons with $E=10^{17}$ eV and $\theta=45^\circ$. The showers were selected so that their $X_{\rm max}$ values are all within less than $1\,\mathrm{g\,/cm^2}$ of $X_{\rm max}=642.5 \, \mathrm{g/cm^2}$. We used a horizontal magnetic field configuration with $|\vec{B}|=50\,\mu$T, parallel to the ground and pointing towards magnetic north. The dashed lines correspond to the average $\langle E_{y}\rangle$ while the shaded area represents the region between $\mathrm{max}(E_{y})$ and $\mathrm{min}(E_{y})$ for ZHAireS (red) and CoREAS (blue) over the 20 simulations. Left panel: $E_y$ component at 50 MHz as a function of the off-axis angle $\alpha$ with respect to the Cherenkov angle $\alpha_C$. The observers are located along the east-west (EW) direction (see Fig. \ref{['fig:esquema']}). Positive coordinates correspond to antenna positions located east of the core. The showers develop from north to south. Right panel: $E_y$ for observers located at distances to the east of the shower core on the ground, corresponding to off-axis angles $\alpha/\alpha_C\simeq 1$ and $\alpha/\alpha_C\simeq 1.5$. The bottom panel of each plot shows the relative difference between the averages, $(\rm ZHAireS - \rm CoREAS)/\rm ZHAireS$, (solid line). The shaded band represents the $1\,\sigma$ spread propagated to the relative difference.
• Figure 4: Same as Fig. \ref{['fig:proton45LDF1']} for primary iron nuclei with $E=10^{17}$ eV and $\theta=75^\circ$ and a vertical magnetic field configuration with $|\vec{B}|=50\,\mu$T, perpendicular to the ground. The simulated showers were selected so that their values of $X_{\rm max}$ are within less than $1\,\mathrm{g/cm^2}$ of $X_{\rm max}=574.5 \, \mathrm{g/cm^2}$. In the left panel the observers are located along the north-south (NS) direction (see Fig. \ref{['fig:esquema']}), with north corresponding to the positive x-axis. In the right panel the antennas are placed at positions north of the shower core.
• Figure 5: Same as Fig. \ref{['fig:proton45LDF1']} for primary electrons with $E=10^{17}$ eV and $\theta=75^\circ$ and a vertical field configuration with $|\vec{B}|=50\,\mu$T, perpendicular to the ground. The 20 simulated showers with each CoREAS and ZHAireS were selected so that their values of $X_{\rm max}$ are within less than $1\,\mathrm{g/cm^2}$ of $X_{\rm max}=715.2 \, \mathrm{g/cm^2}$.
• Figure 6: Peak value of the $y$-component of the electric field ($E_{y,\mathrm{max}}$) at a frequency of 50 MHz as a function of the fraction of the primary energy ($E=10^{17}$ eV) in electromagnetic particles (see text for definition) for protons with $\theta=45^\circ$ (circles), iron $\theta=75^\circ$ (triangles), and electrons $\theta=75^\circ$ (inverted triangles). Filled blue symbols: CoREAS simulations; empty red symbols: ZHAireS simulations.