Reconstructing the cosmic-ray energy from the radio signal measured in one single station

TL;DR

This paper introduces a novel method to reconstruct cosmic-ray energy using radio signals from a single antenna station, achieving better than 15% resolution, which enhances the capabilities of sparse and hybrid radio arrays.

Contribution

A new technique for estimating cosmic-ray energy from a single station's radio data, improving energy resolution and aiding neutrino and cosmic-ray detection.

Findings

Achieves <15% energy resolution with single-station radio measurements.

Utilizes a broad 80-300 MHz frequency bandwidth.

Applicable to hybrid arrays with independent geometry measurements.

Abstract

Short radio pulses can be measured from showers of both high-energy cosmic rays and neutrinos. While commonly several antenna stations are needed to reconstruct the energy of an air shower, we describe a novel method that relies on the radio signal measured in one antenna station only. Exploiting a broad frequency bandwidth of $80-300$ MHz, we obtain a statistical energy resolution of better than 15\% on a realistic Monte Carlo set. This method is both a step towards energy reconstruction from the radio signal of neutrino induced showers, as well as a promising tool for cosmic-ray radio arrays. Especially for hybrid arrays where the air shower geometry is provided by an independent detector, this method provides a precise handle on the energy of the shower even with a sparse array.

Tables

Table 1. Table 1 : Zenith dependence of the parameters A 𝐴 A and s 𝑠 s from Eq. ( 3.2 ) for a station at Moore’s Bay, the South Pole and the site of the Pierre Auger Observatory.

	$A=\left(a\cdot {\theta }^2+b\cdot \theta +c\right)\cdot \frac{{10}^{18}}{m\sqrt{eV}}$			$s=m\cdot \theta +n$
	a	b	c	m	n
Moore’s Bay
	442.46	-281.75	324.96	-0.1584	-0.07943
$\phi >{\phi }_{Cherenkov}$	394.08	-308.36	436.30	0.8070	-1.4098
South Pole
	976.30	-1213.43	626.98	-0.2273	0.05627
$\phi >{\phi }_{Cherenkov}$	643.39	-667.08	478.06	1.3372	-2.1653
Auger
	229.96	-123.75	110.51	-0.1445	-0.09820
$\phi >{\phi }_{Cherenkov}$	214.46	-111.01	119.18	0.5936	-1.1763

Equations

$$\rho (\Delta n)=\frac{{\sum }_i(V_1{)}_i\cdot (V_2{)}_{i-\Delta n}}{\sqrt{{\sum }_i(V_1{)}_i^2}\cdot \sqrt{{\sum }_i(V_2{)}_i^2}}$$

$$\left(\begin{array}{c}{\mathcal{V}}_1(f)\\ {\mathcal{V}}_2(f)\\ ...\\ {\mathcal{V}}_n(f)\end{array}\right)=\left(\begin{array}{cc}{\mathcal{H}}_1^{\theta }(f)& {\mathcal{H}}_1^{\varphi }(f)\\ {\mathcal{H}}_2^{\theta }(f)& {\mathcal{H}}_2^{\varphi }(f)\\ ...& \\ {\mathcal{H}}_n^{\theta }(f)& {\mathcal{H}}_n^{\varphi }(f)\end{array}\right)\left(\begin{array}{c}{\mathcal{E}}^{\theta }(f)\\ {\mathcal{E}}^{\varphi }(f)\end{array}\right)+\text{noise},$$

$$\left(\begin{array}{c}{\mathcal{E}}_{\theta }\\ {\mathcal{E}}_{\varphi }\end{array}\right)=\left(\begin{array}{c}{A}_{\theta }\\ A_{\varphi }\end{array}\right)10^{f\cdot m_f}\exp (\Delta j)$$

$${\Phi }_E={ϵ}_0\cdot c\cdot \Delta t\cdot \sum _i|\vec{E}(t_i){|}^2$$

$$\Phi_{E}^{\prime}=\bigg{(}\frac{1}{\sin^{2}(\alpha)}\Phi_{E,\vec{v}\times\vec{B}}+\Phi_{E,\vec{v}\times(\vec{v}\times\vec{B})}\bigg{)}\cdot\bigg{(}\frac{d_{X_{max}}}{$1\text{\,}\mathrm{k}\mathrm{m}$}\bigg{)}^{2}$$

$$\frac{\sqrt{{\Phi }_E^{\prime }}}{E_{shower}}=A\cdot \exp (-s\cdot (|m_f|\cdot \mathrm{GHz}{)}^{0.8})$$

$$\left(\begin{array}{c}{\mathcal{E}}_{\theta }\\ {\mathcal{E}}_{\varphi }\end{array}\right)=\left(\begin{array}{c}{A}_{\theta }\\ A_{\varphi }\end{array}\right)10^{f\cdot m_f+(f-80\mathrm{MHz}{)}^2\cdot n_f}\exp (\Delta j)$$

$$g^{\prime }=p\cdot g$$

$$g^{\prime}(f)=g(f)\cdot\bigg{(}1+q\cdot\frac{f-\mathrm{200MHz}}{\mathrm{200MHz}}\bigg{)}$$

Full text

1 Introduction and scientific motivation

There is still no commonly accepted answer to the question of the origin of ultra-high energy cosmic rays. Detecting the radio emission of showers may bring us closer to an answer. Either through the detection of air showers, exploiting the excellent energy and composition resolution that comes with the radio method [1, 2, 3, 4, 5], or through the detection of the neutrino counterpart of ultra-high energy cosmic rays. Neutrinos can either be generated in interactions at the sources (astrophysical neutrinos), allowing for multi-messenger detection of for example gamma-rays and neutrinos [6, 7, 8] or as cosmogenic neutrinos [9, 10, 11]. These neutrinos are created through the interaction of ultra-high energy cosmic rays during propagation with the cosmic microwave or other photon backgrounds, making the neutrino flux sensitive to the composition of cosmic rays [12, 13, 14]. All models predict a low flux of these neutrinos beyond ${10}^{16}\text{\,}\mathrm{e}\mathrm{V}$ meaning that current detectors like IceCube are too small to detect a significant flux [15] at the highest energies. Due to the long attenuation length of radio waves of $\mathcal{O}($1\text{,}\mathrm{k}\mathrm{m}$)$ in ice [16, 17, 18] allowing for sparse instrumentation, radio neutrino detectors are a promising alternative.

1.1 Method

When reconstructing particle showers, be it neutrino or cosmic-ray induced, typically two aspects are relevant: the energy and the shower development, i.e. the shower maximum for air showers and the vertex position for neutrinos. Radio detection of air showers has reached maturity already, where the energy is typically derived from reconstructing the footprint of the shower on ground [2, 19, 20, 21]. To this end, the footprint is sampled by many antennas. The radio footprint is typically smaller than the particle footprint, since it is primarily governed by the width of the Cherenkov cone, which depends on the distance to the shower maximum and not the energy of the shower. Especially for events with small zenith angles, this results in measurable signals on only a radius of the order of $100\text{\,}\mathrm{m}$ around the shower axis. In order to sample this footprint with at least three antennas, the radio arrays have to be relatively densely spaced. A much wider spacing would be needed to cost-effectively cover large areas, as it is the case for e.g. the radio component of the proposed AugerPrime upgrade to the Pierre Auger Observatory [22]. At a spacing on the order of $1\text{\,}\mathrm{k}\mathrm{m}$ the chances of being able to sample the full footprint of an air shower are small, unless targeting very horizontal showers ($>60^{\circ}$) [23]. A similar challenge presents itself to proposed radio-based in-ice neutrino detectors (e.g. [24, 25]). With a detector spacing of more than $1\text{\,}\mathrm{k}\mathrm{m}$ it is unlikely for the radio signal from a neutrino interaction to be detected by more than one station. Since radio-neutrino detectors with surface antennas are also sensitive to air showers, detecting cosmic rays and reconstructing their energy will be a valuable tool for validating neutrino arrays, much like using muons in optical-Cherenkov neutrino experiments. Air showers show signal characteristics very similar to neutrino showers, making them both a background and a calibration tool.

In order to fully exploit this calibration signal, radio detection of air showers will need new reconstruction methods for neutrino arrays. In turn, sparse air-shower arrays may also benefit from these methods. We will present a method that derives the energy of an air shower from the signal in a single detector station, exploiting the shape of the frequency spectrum. This method, while carrying merit in itself, can also be seen as step towards energy reconstruction for neutrinos. For the neutrino reconstruction additional complexities such as a much wider Cherenkov angle, frequency-dependent attenuation [16, 17, 18], and a strongly changing index of refraction leading to curved paths will have to be treated. However, we anticipate that using the shape of the frequency spectrum will also be essential for the energy reconstruction of neutrinos.

1.2 Experimental context

Since all large radio neutrino arrays are still in the proposal stage [26, 27, 25, 24], we have developed a method based on CoREAS Monte Carlo simulations [28]. However, in order to embed this method in realistic conditions, we have adopted a station design used in the pilot-stage array of ARIANNA [29], which has already been used for air shower detection [30]. A similar surface component employing log-periodic dipole antennas (LPDAs) is foreseen for all proposed in-ice neutrino arrays, which makes this a solid choice for air-shower reconstruction [24, 25]. We note that the detailed modelling of a detector is not critical for the success of the method, as long as a wide frequency range ($>100$ MHz) is detected.

For the purpose of this article, we adopt a station layout of four upward-facing LPDAs buried in the snow at a depth of roughly $2\text{\,}\mathrm{m}$ below the surface. The antennas are distributed around a central tower, which in practice might be used for power and communications, but simply serves as symmetry center for this example. LPDAs on opposite sides of the central tower are $8\text{\,}\mathrm{m}$ apart from each other and form a pair, as they have the same orientation and therefore measure the same polarization. One pair measures the north-south, the other the east-west polarization of the radio signal.

We adopt typical ARIANNA hardware characteristics. They include a very broad-band LPDA, which is filtered with a $80\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ high-pass and a $500\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ low-pass filter. After amplification the signal is digitized by an analog-digital converter (ADC) with a sampling rate of $1\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$. We also use the ARIANNA-style dual-threshold trigger scheme with a $4\sigma$ coincidence between at least two channels required for the station to trigger. As noise is an important factor for radio arrays, we use recorded noise from the set-up in Moore’s Bay [30]. Adopting ARIANNA characteristics allows us to study the influence of real hardware characteristics as well as background effects to avoid an idealized study of the reconstruction capabilities. In App. A we describe cross-checks for other locations on Earth, to exclude that the height and geomagnetic field at the ARIANNA site are a relevant factor.

In this analysis, we only use data from a single station and neglect possible additional information from other stations detecting the same air shower. While there will be some station coincidences in an array with a spacing in the order of $1\text{\,}\mathrm{k}\mathrm{m}$, especially for horizontal showers, we use this conservative estimate to develop the reconstruction algorithm. The reconstruction of an individual shower will of course improve, if information from two or more stations is provided.

2 Radio Signal Reconstruction

The first step in the event reconstruction is the retrieval of the electric field from the voltage traces measured in the antennas, which in turn have to be retrieved from the recorded ADC counts. To do so, we use the NuRadioReco software framework. A detailed description of NuRadioReco and its event reconstruction methods is presented elsewhere [31], so only the reconstruction of the signal direction and the electric field will be covered here briefly, as they are essential for the energy reconstruction.

2.1 Arrival direction

The distance of $8\text{\,}\mathrm{m}$ between two antennas with the same orientation is small compared to the size of the air shower footprint at ground level. Therefore the radio signal measured in the two channels of an antenna pair can be expected to be very similar, only shifted by the difference in signal travel time between the antennas.

We define the correlation

[TABLE]

between the voltage time traces $V_{1}$ and $V_{2}$ of the two channels in each pair, shifted by $\Delta n$ samples relative to each other. The correlation $\rho(\Delta n)$ is normalized to $-1<\rho(\Delta n)<1$ to make it independent of the signal amplitudes.

Assuming a plane-wave front, we calculate the expected time difference and the corresponding $\Delta n$ between antennas of each channel pair for a given signal direction and calculate $\rho(\Delta n)$. The direction of the recorded signal is determined by finding the direction for which the sum of $\rho(\Delta n)$ of both channel pairs is largest. Using this method, an angular resolution better than $1\text{\,}{}^{\circ}$ is achieved for signals that exceed a signal-to-noise ratio of 4222We define the signal-to-noise ratio (SNR) as the ratio between half the peak-to-peak amplitude and the noise RMS..

2.2 Electric field

With the incoming direction of the radio signal known, it is possible to reconstruct the 3-dimensional electric field trace from the voltage traces measured by the antennas. Since the distances between a station’s antennas are small, we can assume that all four channels measured the same electric field. Then, in Fourier space, the voltage traces in the antennas for a given electric field $\mathcal{E}^{\theta,\phi}$ are given by

[TABLE]

where $\mathcal{H}_{i}^{\theta}$ and $\mathcal{H}_{i}^{\phi}$ are the antenna responses of the ith antenna to the $\theta$ and $\phi$ components of the electric field. Eq. (2.2) is an overdetermined system of linear equations that can be solved for $\mathcal{E}^{\theta}(f)$ and $\mathcal{E}^{\phi}(f)$ via a $\chi^{2}$ optimization to reconstruct the electric field. However, this method does not respond well to noise in the voltage traces [31]. We therefore use the so-called forward folding method instead, which provides a more accurate electric field reconstruction, especially with regards to the frequency spectrum by using an analytic model of the radio signal.

In the frequency domain, the electric field can be approximated by the equation

[TABLE]

that depends on only four parameters: The signal amplitudes $A_{\theta}$, $A_{\phi}$ of the two orthogonal polarizations, the frequency slope $m_{f}$ and a frequency-independent phase offset $\Delta$. $f$ is the frequency and $j$ the imaginary unit.

The electric field is reconstructed by calculating $\mathcal{V}_{i}(f)$ for a given $\mathcal{E}_{\theta}(f)$, $\mathcal{E}_{\phi}(f)$ using Eq. (2.2) and optimizing $A_{\theta}$, $A_{\phi}$, $m_{f}$ and $\Delta$ through a $\chi^{2}$ fit to the recorded voltage time traces. An example of such a forward folding reconstruction is displayed in Fig. 1. The voltage traces were generated by folding the electric field from an air shower simulated using CoREAS [28] with the detector response and adding a noise trace recorded by an ARIANNA station in the field. It shows that the forward folding method avoids the problem of overestimating of the electric field at higher frequencies that the standard method suffers from, resulting in a more accurate spectrum reconstruction. Especially the reconstruction of the frequency slope $m_{f}$ is much more accurate than with the standard method. For a detailed comparison between the performance of the two methods, see [31].

From the electric field trace we calculate the energy in the radio signal per unit area, called the energy fluence [2]

[TABLE]

where $\epsilon_{0}$ is the vacuum permittivity, $c$ is the vacuum speed of light and $\Delta t$ is the sampling rate. We note that we do not need to subtract noise from the energy fluence, because the forward-folding technique reconstructs an electric field free of noise.

3 Energy estimator

Our goal is to use the energy fluence $\Phi_{E}$ as an estimator in order to measure the energy in the electromagnetic shower. In our scenario, each air shower is only detected by a single radio station and no additional information about the shower geometry is available. Thus, an energy reconstruction from the energy fluence at a single point is challenging, as it is influenced by a number of factors other than the shower energy:

Angle to the magnetic field:

The dominant emission process for air showers is the so-called geomagnetic emission, whose energy fluence is proportional to $\sin^{2}(\alpha)$ where $\alpha$ is the angle between the shower axis and the geomagnetic field.

Distance to the shower maximum:

The radio signal gets weaker with increasing distance between the shower and the radio antenna. Since the largest part of the radio signal is emitted near the shower maximum, the signal energy fluence roughly scales with the geometric distance to the shower maximum $d_{X_{max}}$ squared.

Viewing Angle:

We define the viewing angle $\varphi$ as the angle between the shower axis and a line from the shower maximum to the antenna. Emission is strongest if $\varphi$ is close to the Cherenkov angle (defined by $\cos(\varphi_{Cherenkov})=\frac{1}{n}$, where $n$ is the index of refraction at the shower maximum) and becomes weaker the larger the difference between the two gets.

For an estimator to achieve a good energy resolution, we need to take these factors into account and correct for them.

3.1 Procedure

Correcting for the angle to the magnetic field is rather straightforward, since the shower axis is very close to the signal direction and the magnetic field at any site can be measured directly. For the correction, it is useful to introduce a coordinate system that has one axis aligned with the shower axis $\vec{v}$, one in the direction of the Lorentz force $\vec{v}\times\vec{B}$ and the third one perpendicular to those two $\vec{v}\times(\vec{v}\times\vec{B})$. The energy fluence can then be split into the energy fluence of the two polarizations $\Phi_{\vec{v}\times\vec{B}}$ and $\Phi_{\vec{v}\times(\vec{v}\times\vec{B})}$. Since only the geomagnetic emission, which is polarized in the $\vec{v}\times\vec{B}$ direction, is affected by the angle to the magnetic field, only $\Phi_{\vec{v}\times\vec{B}}$ is corrected for this effect.

The distance $d_{X_{max}}$ to the shower maximum is harder to account for, because the atmospheric depth $X_{max}$ of the shower maximum is unknown. Fortunately, however, this effect is most relevant for inclined showers, where the shower maximum is far away from the observer. In this case, event-by-event fluctuations of $X_{max}$ have little effect and $d_{X_{max}}$ is mainly a function of the zenith angle [32]. It is therefore sufficient to assume a typical value for $X_{max}$ and calculate the expected $d_{X_{max}}$ based on the shower inclination. Based on [33] we use an average $X_{max}$ of $750\text{\,}\frac{g}{cm^{2}}$333In fact, the value chosen for $X_{max}$ has very little effect on the energy reconstruction. If a different $X_{max}$ had been chosen, the zenith-dependent parameterization of the effect of the viewing angle (Sec. 3.2) would compensate this.. To make sure that this assumptions was correct, we checked for correlations between the uncertainties on the reconstructed energy and $X_{max}$ and found no such biases.

This lets us correct for those two effects by defining a corrected energy fluence

[TABLE]

which will be used as the basis for the energy estimator. We note that we ignore the influence of the second order charge-excess emission process on the $\Phi_{E,\vec{v}\times\vec{B}}$-component as it is small enough to not limit the resulting energy resolution. Since geomagnetic and charge excess emission can interfere both constructively and destructively depending on where the observer is on the $\vec{v}\times\vec{B}$ axis, the effect of the charge excess emission will largely average out and not lead to a systematic offset. In the future, the term $\sin^{2}(\alpha)$ may be replaced by $a^{2}+(1-a^{2})\sin^{2}(\alpha)$ following the same approach as in [34] where $a$ is the relative strength of the charge-excess component that depends on the distance to $X_{max}$ and the shower axis, and can be parameterized as function of zenith angle using the work of [32].

The third correction for the effect of the viewing angle is more complex and has to be assessed in a Monte Carlo study. We use a first set of simulated air showers to derive the correction factor for the viewing angle and a second set to test the full energy reconstruction (see Sec. 4).

The Monte Carlo data set was produced using 545 CoREAS [28] simulations of proton-induced air showers and their radio emission. The showers cover an energy range of ${10}^{16}\text{\,}\mathrm{-}$ $5\text{\times}{10}^{19}\text{\,}\mathrm{e}\mathrm{V}$ and random arrival with zenith angles between $0\text{\,}{}^{\circ}$ and $80\text{\,}{}^{\circ}$ and a uniform azimuth distribution. The energy range is set by the typical threshold of radio signals above the thermal and Galactic noise floor. This depends in detail on the experimental set-up but is typically at several PeV.

For each shower, the radio signal was simulated at 160 positions arranged in a star-like pattern around the shower core at an altitude of $30\text{\,}\mathrm{m}$ above sea level (the height of the ARIANNA detector on the Ross Ice shelf), and the detector response to the radio signal at each position was simulated. Radio simulations are in general time consuming, but every observer position can be used independently providing additional statistics. Also, since the radio signal stems from the electromagnetic component of the air shower, which is in general rather smooth, small air shower sets provide enough information. We find that we are not statistics limited (see also App. A).

A simple-threshold trigger was simulated by requiring that at least two channels recorded a pulse larger than $60\text{\,}\mathrm{m}\mathrm{V}$ ($\sim 3V_{rms}$) for the reconstruction to continue. A 10th order Butterworth bandpass filter with a passband of 80\text{\,}\mathrm{-}$$300\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} was applied and a full event reconstruction was performed, using the forward folding technique described in Sec. 1. The resulting data set was used to develop an energy estimator that can also account for the effect of the viewing angle.

3.2 Correction for the viewing angle

The shape and amplitude of a radio pulse emitted by an air shower differs depending on the viewing angle $\varphi$ under which it is seen. Under small viewing angles, signals from different stages in the shower’s development arrive simultaneously, leading to an amplified emission over a wide frequency band [35]. The resulting dependency of the frequency spectrum on $\varphi$ has already been observed in the $30-80$ MHz band [36, 37, 38, 23] and is used to reconstruct air shower energies by the ANITA experiment [39]. The method used by ANITA, however, requires dedicated Monte Carlo simulations for each reconstructed event and is therefore not suitable for large arrays with regular detections.

The radio signal is strongest if the viewing angle coincides with the Cherenkov angle at the location of the shower maximum. In this case, the radio signals emitted at the shower maximum add up coherently, leading to an amplification of the signal. If the event is seen from further off the Cherenkov cone, overlapping signals lose coherence and the amplitude decreases, as can be seen in the top row of Fig. 2. However, shorter wavelengths lose coherence faster than longer ones. The effect this has on the frequency spectrum can be observed with the frequency slope parameter $m_{f}$ (see Eq. 2.3), shown in the bottom row of Fig. 2. The small absolute value of $m_{f}$ at the Cherenkov angle corresponds to a relatively flat frequency spectrum, while the larger $m_{f}$ off the Cherenkov cone corresponds to a sharp dropoff towards higher frequencies in the spectrum. The event-by-event fluctuations in $\sqrt{\Phi_{E}^{\prime}}/E_{shower}$ and $m_{f}$ are large for $\theta<10^{\circ}$, but they decrease with increasing zenith angle. It is also worth noting that neither the change in $\sqrt{\Phi_{E}^{\prime}}$ nor in $m_{f}$ is symmetric around the Cherenkov angle.

Since $m_{f}$ is a function of $\varphi-\varphi_{Cherenkov}$, it is possible to use it as a proxy for the viewing angle. This is done in Fig. 3, which shows $\sqrt{\Phi_{E}^{\prime}}/E_{shower}$ as a function of $m_{f}$, with small $|m_{f}|$ corresponding to events seen by observers close to the Cherenkov ring and larger $|m_{f}|$ to those seen from further off the Cherenkov angle. Events are split into two groups: Those where the station is located inside the Cherenkov ring and those seen from outside of it. While $\sqrt{\Phi_{E}^{\prime}}/E_{shower}$ and $m_{f}$ show a clear correlation for larger zenith angles, the scatter for the $\theta<10^{\circ}$ events is very large. This is likely due to the small angle between the shower axis and the geomagnetic field, which causes the charge-excess emission to become more prevalent compared to the geomagnetic emission. Since both signal components have a different polarization, they interfere constructively at positions in positive $\vec{v}\times\vec{B}$ direction and destructively in negative $\vec{v}\times\vec{B}$ directions. This causes the lateral distribution function (LDF) of the radio signal to no longer be rotationally symmetric in the shower plane. Together with the smaller dependence on the viewing angle (Fig. 2), this leads to the large scatter (see also discussion in App. A).

We divide the events into 20 zenith angle bins with equal sky coverage and fit the function

[TABLE]

to the distribution in each bin separately for events seen from within and from outside of the Cherenkov cone. The results are shown in Fig. 4. The zenith dependence of $A$ can be parameterized by a second-order polynomial, the one of $s$ by a straight line. The specific values for this parametrization are shown in Tab. 1.

Using these parameterizations, Eq. (3.2) can now be solved for $E_{shower}$ to obtain an energy estimator using the zenith angle $\theta$ and the radio energy fluence $\Phi_{E}$ and frequency slope $m_{f}$. However, since events viewed from inside the Cherenkov cone have to be treated separately from those seen from outside of it, a way to distinguish between the two is needed, which will be addressed in the next section.

3.3 Inside vs. outside of the Cherenkov cone

When reconstructing the electric field in Sec. 1, we approximated frequency dependence of $\log_{10}(|\vec{\mathcal{E}}(f)|)$ by a straight line. As is shown in Fig. 5, the actual radio signal deviates from this, with the frequency spectrum curving downward at higher frequencies for events seen from inside the Cherenkov ring and upward for those seen from outside of it. While the difference is small enough that it does not need to be taken into account when reconstructing $\Phi_{E}$ and $m_{f}$, it can be used to distinguish between the two event categories.

We extend the reconstruction of the electric field by an additional step, in which we modify Eq. (2.3) by introducing an additional 2nd order term in the exponent:

[TABLE]

This pulse is fitted to the recorded voltage traces as described in Sec. 1, using the results from fitting Eq. (2.3) as starting values. The amplitudes $A_{\theta}$ and $A_{\phi}$ are allowed to change, while $m_{f}$ is kept constant.

The relation between $\varphi$, $m_{f}$ and $n_{f}$ is shown in Fig. 6. Showers seen from within the Cherenkov ring tend to have a smaller $n_{f}$ than those seen from outside of it. A simple linear function (black dotted line) can be used to separate the events into two groups: Those above the line and those below it. Several cuts have been tested leading to similar efficiency and purity, so we have chosen to use the simplest. The histogram on the right shows the viewing angles relative to the Cherenkov cone of the two groups, demonstrating that these provide a good separation between events viewed from within and outside of the Cherenkov ring. A fraction of 90% of the events are either identified correctly or are less than $0.05\text{\,}{}^{\circ}$ away from the Cherenkov angle. A mistake in the separation of events this close to the Cherenkov ring will result only in minor uncertainties on the reconstructed energy, as the energy parametrizations for both cases are very similar near the Cherenkov angle. One may notice that for large $|m_{f}|$, the discriminator seems to work very poorly, as events outside the Cherenkov ring tend to have small values of $n_{f}$. In principle, this could be improved by also defining a maximum value for $|m_{f}|$ above which all events are classified as being outside the Cherenkov ring. In practice however, the signals from these events are so weak that they are unlikely to be detected, making this additional step unnecessary.

As it was the case for the energy parametrization, this method works best for inclined showers. The difference in $n_{f}$ for showers seen from within and outside of the Cherenkov ring tends to increase with the zenith angle. For small zenith angles, both event categories mostly overlap.

Knowing if the station is inside or outside of the Cherenkov ring, Eq. (3.2) can now be used to determine the electromagnetic shower energy. The resulting energy reconstruction algorithm has been implemented as the cosmicRayEnergyReconstructor module in the NuRadioReco software framework [31], where one can see the full implementation in detail and use it to reconstruct the energy accordingly [40].

4 Obtainable energy resolution

Having developed an energy parametrization and a way to distinguish between air showers seen from inside and outside the Cherenkov ring, these can now be combined to reconstruct the shower energy. By performing this energy reconstruction on simulated air showers and comparing the result to the true energy, we can determine its performance under realistic conditions.

4.1 Monte Carlo data-set

To estimate the uncertainty of the reconstructed energies, a set of 200 CoREAS [28] simulations of air showers and their radio emissions was used. To mimic a realistic data set, the showers had an isotropic distribution of incoming directions up to a zenith angle of $80^{\circ}$ and energies in the range of {10}^{17}\text{\,}\mathrm{e}\mathrm{V}$<E<${10}^{20}\text{\,}\mathrm{e}\mathrm{V} with a spectral index of $-2$. For each air shower, the radio signal was simulated at 160 positions, arranged in a star-like pattern around the shower core. To achieve a realistic distribution of station positions, for each simulated air shower, 20 random positions within a radius of $2\text{\,}\mathrm{k}\mathrm{m}$ around the shower core were produced. For each of these random positions, the simulated electric field closest to it was determined and used in the analysis. More detailed information on the used Monte Carlo set is provided in App. B.

Simulated electric field traces were folded with the antenna response and the effect of the signal chain was simulated to obtain the voltage traces. To account for the effect of noise, one of 100,000 forced-trigger events recorded by ARIANNA stations were randomly selected and its recorded voltage trace added. A dual-threshold trigger [29] was simulated by requiring the voltage pass both an upper and a lower threshold of $\pm$80\text{,}\mathrm{m}\mathrm{V} (about 4 times the typical noise RMS) within $5\text{\,}\mathrm{n}\mathrm{s}$ in at least 2 channels. Of the 4,000 radio pulses, 597 passed this selection. A 10th order Butterworth filter with a pass-band of $80\text{\,}\mathrm{-}300\text{,}\mathrm{M}\mathrm{H}\mathrm{z}$ was applied to the voltages and a reconstruction of the incoming direction and the electric field was performed. The electromagnetic shower energy was estimated from Eq. (3.2). Since the electric field reconstruction corrects for the detector response to the radio signal and the same detector description was used for simulation and reconstruction, the effect of uncertainties on the hardware components is neglected here. These concerns are addresses separately in Sec. 4.3.

4.2 Energy resolution

The resulting energy resolution is shown in Fig. 7. While the energy reconstruction does not work well for vertical showers, for more inclined showers the distribution has a median of $-0.01^{+0.21}_{-0.11}$ with the uncertainties representing the 68% quantile. While the uncertainty is small for most events, there are several outliers where the energy was greatly overestimated.

Most of these outliers happen with events that have a low signal-to-noise ratio. Because some periods have a higher noise level than normal [41], even a trigger with a threshold of $80\text{\,}\mathrm{m}\mathrm{V}$ can produce events with a small signal-to-noise ration. Another problem can occur with the direction reconstruction. Because each channel pair measures a different polarization, it is possible for the radio pulse to only be detected in one pair. In this case, the direction reconstruction is likely to produce a wrong result which in turn leads to a wrong reconstruction of the electric field and the energy. We therefore impose a stricter cut requiring that at least one channel in each pair has as $\mathrm{SNR}>4$. The resulting energy resolution is shown in Fig. 8 (left). With this cut, most of the outliers disappear and the resolution improves, with a median of $-0.02^{+0.12}_{-0.10}$.

It turns out that most of the outliers in Fig. 7 are caused by events viewed from within the Cherenkov ring being mistaken for events seen from outside of it. While an ARIANNA-like detector has no other way to identify if the station was inside or outside the Cherenkov ring, this task may be trivial for other radio detector designs, for example those that are part of a hybrid detector. Fig. 8 shows the energy resolution if reliable information if the shower was seen from inside the Cherenkov ring or outside of it was available from some other source. In this case most outliers disappear and the resulting distribution has a median of $-0.01^{+0.15}_{-0.11}$ with the uncertainties representing the 68% quantile, having applied no SNR-cut. The remaining outliers are still mostly events with a relatively low signal-to-noise ratio, in which the electric field reconstruction did not converge properly.

It should be noted that this method, like all radio methods, reconstructs the energy in the electromagnetic cascade of the shower [34]. To get the cosmic ray energy, the invisible energy, which can be up to $\sim 20\%$ of the total shower energy, has to be taken into account as well [42].

4.3 Systematic effects of hardware uncertainties

Up to this point, the event simulation and reconstruction have used the same detector description, implying perfect knowledge of the detector response. This is, of course, not realistic, so in this section we will investigate the effect of hardware uncertainties on the energy resolution.

The procedure is as follows: We repeat the same simulation and reconstruction as in Sec. 4, but when the detector response is simulated, certain aspects of the detector description are modified. Then, we reconstruct the event using the original detector description. This mimics the effect that uncertainties in the hardware would have on the event reconstruction.

The first investigated uncertainty is a deviation in the orientation of the antennas. We produce two sets of simulated events with one of the antennas rotated by $2^{\circ}$ and $5^{\circ}$ around the z-axis and perform a full event reconstruction using the default antenna orientation on both sets. The uncertainties on the reconstructed energies show no significant difference to those in Fig. 7. We therefore conclude that the influence of uncertainties on the antenna orientation on the reconstructed energy is negligible since $5^{\circ}$ precision in the antenna orientation is easily achievable.

For the effect of uncertainties on the antenna or amplifier gain $g(f)$, we consider two different effects: The gain may be off by a constant factor for all frequencies. In that case, the energy fluence of the reconstructed radio signal would change, but the frequency slope $m_{f}$ would remain unchanged. To test this, we scale the gain $g$ by a constant factor $p$:

[TABLE]

If the uncertainties in the amplifier gain are frequency-dependent, they may influence the reconstructed frequency slope $m_{f}$. To investigate this, the gain is modified by

[TABLE]

Because we operate in the 80\text{\,}\mathrm{-}$$300\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} frequency band, this means the gain will be increased at one edge of the frequency band and decreased at the other by the same factor $q/2$. This modification will mainly affect the reconstruction of $m_{f}$, while the effect on $\Phi_{E}$ should be small.

The effect of these modifications to the amplifier gain are shown in Fig. 9. On the left, the gain was scaled by $p=0.9$, on the right the slope was modified by $q=0.1$. In both cases, the median shifted by $\approx 10\%$ compared Fig. 7.

We tested changes in the amplifier gain for several values of p in the range $0.9<p<1.1$ and q in the range $-0.1<q<0.1$. In both cases, the shift of the median of the energy resolution is roughly linear with respect to p and q. The width of the Gaussian function fitted to the histogram showed little change, but the upper limit of the 68% quantiles increased for some modified amplifier gains. This was caused by errors in determining if the air shower was seen from inside or outside the Cherenkov ring, which caused additional outliers. Under the assumption of a known shower geometry, the 68% quantiles remained constant.

5 Discussion

We demonstrated that the signal detected by a single radio station is sufficient to accurately measure the energy in the electromagnetic part of an air shower. While the statistical uncertainty is better than 15%, additional systematic uncertainties from the hardware response, namely their influence on the energy fluence and the frequency spectrum, have to be considered, underlining the importance of a thorough calibration.

The estimator was developed in the context of a surface station of a neutrino array, but the concept can be applied to different detectors, both for cosmic rays and neutrinos, as long as the signal is measured in a broad-band going beyond $100\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$.

5.1 Implications for cosmic-ray detectors

Our results may be applied to dedicated radio-based cosmic ray detectors in order to improve the accuracy of energy measurements. However, it requires a larger frequency band than the 30\text{\,}\mathrm{-}$$80\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} band that most radio detectors today operate in, for two reasons: A larger frequency band makes the measurement of the frequency slope parameter $m_{f}$ easier, as simply more data points are available for the reconstruction. Also, in the 30\text{\,}\mathrm{-}$$80\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} band, the frequency spectrum can not be described by a simple exponential function [37, 36], so a more complex parameterization would have to be used, which makes the reconstruction more challenging. Together with an improved signal to noise ratio at higher frequencies [43], this makes the use of frequency bands outside the 30\text{\,}\mathrm{-}$$80\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} band an attractive option.

The method developed in this paper is especially well suited for sparse arrays (antenna spacing $\sim\mathcal{O}$(km)), where many events are only detected by a small number of stations. Sparse arrays also often focus on the detection of inclined air showers, for which this method yields the best results.

In addition, combining the radio measurement of a shower with the measurement of muons provides an excellent estimator for the primary particle. It is shown in [44] that the ratio of the energy in the electromagnetic component of the shower, represented by the radio energy estimator, and the hadronic component, derived from the muonic component can used as powerful discriminator. Thus, building a hybrid detector of radio antennas and muon detector, may successfully address the composition problem. If using the energy reconstruction method presented here, the array can even be sparse and therefore low-cost and competitive. It is also worth noting that one of the major sources for uncertainties is confusion between showers seen from inside and outside the Cherenkov ring. For hybrid experiments combining radio with other detection methods, this issue may be trivial (see e.g. [45]). We also developed and parameterized this method for a detector located at the South Pole and at the site of the Pierre Auger Observatory in Argentina (see App. A).

5.2 Implications for neutrino detectors

The ability to measure the shower energy is a crucial requirement for any radio-based neutrino observatory. Experiments attempting to detect neutrino interactions in ice space their stations so far apart that a signal is expected to be only detected by a single station in most cases. This makes the ability to reconstruct the energy of a particle cascade from the measurement by one station alone essential.

By using air showers as a test case, we have explored a feasible way to achieve this. Because of the similarities between radio signals from air showers to those from neutrino-induced cascades in the ice, it is to be expected that our approach can be applied to neutrino detection as well. However, the radio signal from neutrinos shows additional complexities, such as an unknown vertex position, propagation effects through the ice and no dominant signal polarization due to a negligible geomagnetic emission. Therefore, further development is necessary in order to apply this method to neutrino detections. However, this method illustrates how a large band-width can be instrumental for the reconstruction and that methods like the forward folding of the pulse are robust to noise contributions.

6 Conclusions

We have presented a method to reconstruct the energy in the electromagnetic cascade of an air shower using radio signals recorded in the $80-300$ MHz band by a single detector station. By using a forward folding technique, we are able to measure the frequency spectrum of the signal more accurately and use it to estimate the viewing angle in relation to the Cherenkov angle. This information is combined with the energy fluence and the signal direction to provide an energy estimator.

The accuracy of this method under realistic noise conditions and using a realistic triggering scheme was investigated. For events with a sufficiently high signal-to-noise ratio, statistical uncertainties around 10% are achievable for the energy. Uncertainties in the orientation of the antennas were shown to be negligible while uncertainties in the amplifier gain have a significant effect and propagate linearly as systematic uncertainties on the reconstructed shower energy.

For cosmic-ray detectors, this method can be used to improve the accuracy of energy measurements, especially in sparse arrays, where often only measurements from a small number of stations are available. Regarding neutrino detectors, we have shown that it is possible to reconstruct the energy of a particle cascade from the signal of a single radio station, indicating a possible route towards the energy reconstruction for neutrinos at the highest energies.

7 Acknowledgements

We are thankful for the fruitful collaboration with our colleagues from the ARA and ARIANNA collaborations, especially for providing us with measured noise for this analysis. We acknowledge funding from the German research foundation DFG under grants NE 2031/2-1 and GL 914/1-1.

Appendix A Parameterization for other locations on Earth

Two additional tests were performed for other locations on Earth. The first location was the South Pole, which is the location of IceCube and its future extension IceCube-Gen2, which is planned to include a large radio array. Compared to Moore’s Bay, the South Pole is at a higher elevation of $2800\text{\,}\mathrm{m}$ above sea level. While the geomagnetic field is nearly vertical at both the South Pole and Moore’s Bay, it is weaker at the South Pole. Fig. 10 shows the dependence of the corrected energy fluence of the radio signal $\Phi_{E}^{\prime}$ and spectral slope $m_{f}$ for different zenith angles for a radio detector at the South Pole. The distribution looks similar to the one shown in Fig. 2, indicating that the same method will work. Because of the higher elevation, the shower maximum is much closer to the detector than at Moore’s Bay and in some cases, the shower can even be clipped by the ground before reaching its maximum. Shower clipping is a known issue for the radio detection of air showers [2] and reduces the achievable accuracy, as an unknown fraction of energy has not been deposited in the shower.

The second location we simulated was the site of the Pierre Auger Observatory, where the Auger Engineering Radio Array (AERA) is located [2]. Also, the Auger Observatory is in the process of being upgraded with radio antennas at all surface detectors [22]. The two main differences compared to Moore’s Bay are the higher elevation of roughly $1500\text{\,}\mathrm{m}$ and the geomagnetic field, which has a zenith angle of $54\text{\,}{}^{\circ}$ and is weaker compared to Moore’s Bay. Fig. 11 shows the dependence of the corrected energy fluence $\Phi_{E}^{\prime}$ and spectral slope $m_{f}$ on the viewing angle for a station at an elevation of $1564\text{\,}\mathrm{m}$, where the Pierre Auger Observatory is, and at sea level (same longitude and latitude). As one can see, the differences in elevation only have a minor influence. The overall trend is again the same, one even observes less scatter at small zenith angles.

For a station in Antarctica, the energy of an air shower with a small zenith angle is difficult to reconstruct, because two effects coincide: The steep and featureless lateral distribution at small zenith angles and an almost vertical magnetic field. For small zenith angles, small angles $\alpha$ to the geomagnetic field cause the radio signal to no longer be rotationally symmetric around the shower axis. In addition, the lateral distribution of $\Phi_{E}^{\prime}$ and $m_{f}$ is flatter and provides less discrimination power for the position with respect to the shower axis. Because of the inclined magnetic field at the Auger site, the scatter at small zenith angles is smaller. Small values of $\alpha$ instead occur for more inclined showers, where they pose a smaller problem because of the more pronounced dependence of $\Phi_{E}^{\prime}$ and $m_{f}$ on the viewing angle.

The parametrizations shown in Eq. (3.2) and Fig. 4 have been calculated for the South Pole and the Auger site as well and are shown in Tab. 1 along with the ones for Moore’s Bay.

A detailed investigation of the achievable energy resolution at these locations would go beyond the scope of this paper, but it stands to reason that our method also works elsewhere and might even be able to provide better energy measurements for small zenith angles at locations with a less vertical magnetic field.

Appendix B Monte Carlo Data set

The Monte Carlo data set to evaluate the performance of the energy estimator was generated by randomly generating an incoming direction for the cosmic ray from an isotropic distribution and an energy from a distribution with a spectral index of $-2$. Out of a set of CoREAS files, the one with properties most similar to those generated was selected for each generated event. The resulting distribution of zenith angles, azimuths and energies is shown in Fig. 12.

In each individual CoREAS file, observers are placed on a star-shape like pattern. Observers are placed symmetrically around the shower axis on the $\vec{v}\times\vec{B}-$ and the $\vec{v}\times(\vec{v}\times\vec{B})-$ axis, as well as on two lines bisecting the angle between the two axes.