Constraints on the Diffuse Flux of Ultra-High Energy Neutrinos
from Four Years of Askaryan Radio Array Data in Two Stations

Ultra-high energy neutrinos (UHE, $>10^{17}$ eV) are a unique window on the distant, high energy universe. In addition to gravitational waves, they are the only Standard Model messengers capable of traveling cosmic distances undeflected and unattenuated. Cosmic rays have their trajectories bent by magnetic fields, and for sources more distant than ${\sim}50$ Mpc, above ${\sim}10^{19.5}$ eV cosmic rays are expected to be degraded in energy through interactions with the Cosmic Microwave Background (CMB) via the Greisen-Zatsepin-Kuz’min (GZK) effect [1, 2]. Cosmic ray nuclei are additionally degraded in-flight to earth through their natural beta and inverse-beta decay processes, as well as photo-disintegration e.g. the Giant Dipole Resonance [3]. High-energy gamma rays ($\gtrsim 100$ TeV) are similarly expected to pair-annihilate off the CMB and Extragalactic Background Light (EBL) [4].

Predictions for the sources of very high energy neutrinos fall broadly into two classes. First, astrophysical neutrinos are expected from the site of cosmic ray acceleration, for example gamma ray bursts and active galactic nuclei [5, 6]. The IceCube experiment has confirmed the existence, and measured the spectrum, of TeV-PeV astrophysical neutrinos [7], and has identified a first potential source in the blazar TXS 0506+056 [8, 9]. Second, cosmogenic neutrinos are expected from the destruction of cosmic rays through the aforementioned processes [10]. A more complete discussion of how the flux of cosmogenic neutrinos depends on the primary cosmic-ray composition, and the effects of various interaction and decay processes, can be found in the literature [11, 12, 13, 14, 15].

At energies above above $10^{16}$ eV, low predicted fluxes [16, 17] combined with small expected cross sections [18, 19] lead to $\mathcal{O}(10^{-2})$ neutrino interactions per cubic-kilometer of ice per year per energy decade. As such, the active volumes of the instruments required to detect this UHE flux must necessarily approach the scale of $100$ km³ water equivalent. Several experiments are operating or under-construction to reach this high energy flux, including IceCube [20], Pierre Auger [21], NuMoon [22], ANITA [23], ARIANNA [24], GRAND [25], and ARA [26], which is the focus on this work.

The Askaryan Radio Array (ARA) is a UHE neutrino detector deployed at the South Pole seeking to observe these ultra-high energy neutrinos. ARA searches for neutrinos by looking for the broadband (few hundred MHz to few GHz) radio impulse, or “Askaryan emission” [27, 28], that accompanies neutrino-nucleon interactions. This effect, caused by a $\sim 20\%$ negative charge asymmetry in electromagnetic showers in media, and acting as a coherently radiating current distribution, has been observed in the laboratory at accelerator facilities [29]. The radiation has a Cherenkov-like beam pattern, with a cone thickness of a few degrees. The leading edge of the electric field pulse points toward the shower axis. Experiments looking for Askaryan radiation are deployed in dielectric media such as ice, salt, and sand, which are expected to be sufficiently transparent to radio waves as to make the radio signal observable. In the case of ARA, the long (generally greater than 500 m [30]) attenuation length of radio waves in South Pole ice allows naturally occurring detector volumes to be instrumented sparsely and economically. A diagram of how a neutrino interaction might be observed in an ARA detector is given in Fig. 1.

In this paper, we report constraints on the diffuse flux of UHE neutrinos over the energy interval $10^{16}-10^{21}$ eV. This result is based on two complementary searches for neutrinos in four years of data from ARA stations A2 and A3 recorded between February of 2013 and December of 2016. This paper is organized as follows. In Sec. II, we describe the ARA instrument. In Sec. III, we describe the data analysis methods used in two parallel analyses, and in Sec. IV we discuss our findings. In Sec. V, we discuss systematic uncertainties. Finally, Sec. VI we discuss the result and its implications, as well as prospects for the future. We also include an appendix, App. A, where we discuss the calculation of our limit and detail the livetime of the instrument.

The Askaryan Radio Array is a UHE radio neutrino detector consisting of five stations located a few kilometers grid-west of the geographic South Pole in Antarctica, as drawn in Fig. 2 [26]. A single station consists of 16 antennas, eight for detecting horizontally-polarized (HPol) radiation and eight for detecting vertically-polarized (VPol) radiation, along with signal conditioning and Data Acquisition (DAQ) electronics. The antennas are deployed at the bottom of holes at up to 200 m depth on four “measurement strings,” forming an rectangular solid 20 m tall and 15 m deep and wide. At each corner of the rectangle an HPol quad-slotted cylinder antenna sits a few meters above a VPol wire-frame bicone antenna. Each antenna is approximately sensitive to radiation in the 150-850 MHz band [26]. Two “calibration strings” are deployed about 40 m radially away from the center of the station. Each calibration string contains a VPol and an HPol antenna, and is capable of emitting broadband RF pulses, which provide an in-situ calibration of station geometry and timing, as well as a measurement of livetime.

Construction of ARA began in 2011, when a prototype station (Testbed) was deployed [26, 31] at 30 m depth to evaluate the RF environment and electronics performance. The first design station (A1) was deployed in 2012, but only up to 100 m depth due to limited drill performance. In 2013, two deep stations (A2, A3) that are the focus of this work were deployed at up to 200 m depth [32]. Two more 200 m depth stations (A4, A5) were deployed in 2018.

Our data analysis searches for a diffuse flux of neutrinos between $10^{16}{-}10^{21}$ eV. The analysis is designed to remove background events, principally thermal and anthropogenic noise, while preserving sensitivity to neutrinos. The analysis proceeds in a “blind” fashion, where the ARA data is divided into two subsets. A “burn” sample of 10% of the data, which is assumed to be representative of the full data sample, is set aside and used to tune cuts and understand backgrounds. The remaining 90% of the data is kept blinded until cuts are finalized. Before unblinding, it was decided that in the absence of a detection, the analysis with the best expected limit would be used to set the limit.

After rejecting data containing bursts of surface activity, both analyses examined the neutrino signal region in A2 before examining the signal region in A3. Each analyses’ individual unblinding results are discussed below in Sec. IV.1 and Sec. IV.2.

Neither analysis observes a statistically significant excess of events, observing zero events against the estimated background. In the absence of detection, in Fig. 8 we compute the 90% confidence level (CL) upper limit on the diffuse flux of neutrinos.

Further details on the upper limit calculation, including inclusion of the systematic uncertainties discussed in Sec. V, can be found in App. A. Inclusion of the systematic uncertainties in the limit has an $\mathcal{O}(5\%)$ effect. We report the limit set by Analysis A, which had slightly superior expected sensitivity, by up to 15%, depending on energy. As a benchmark, the number of events expected to be observed in the analysis ranges from 0.25 for an unbroken extrapolation of the astrophysical neutrino flux as measured by IceCube with a spectral index of -2.13 [46], to 0.027 in the case of a cosmogenic neutrino flux where protons make up only 10% of cosmic ray primaries [16].

In Fig. 7, we present the analysis efficiency of Analysis A for both A2 and A3; we plot the average signal efficiency, taking into account the variations due to different run configurations and their respective livetimes. The signal efficiency is calculated by simulating neutrinos in AraSim, and taking the ratio of the number of neutrinos passing the analysis cuts to the number of neutrinos that trigger the detector. We show the efficiency for Analysis A, as it is the analysis used to set our limit, though the efficiencies for Analysis B (which was developed in a parallel and independent fashion) are comparable. On the left, we show the efficiency as a function of SNR, where SNR is computed as the third highest $V_{peak}/RMS$, where $V_{peak}$ is the highest absolute voltage peak in a waveform, and the RMS is the root-mean-square of the voltage values in that waveform. We present the figure with this definition of SNR as it more closely aligns with that commonly used for comparison purposes in the literature. The analysis becomes efficient near an SNR of 6, and does not fully saturate to a value between 75-90% until it is above an SNR of 8. The saturated efficiency for A3 is $\sim 10$% lower than for A2 because A3 required a larger angular cut region to reject surface events, as discuss in Sec. III.4. On the right, we show the efficiency as a function of energy. At $10^{16}$ eV, the analysis has a relatively low efficiency of about 5%. The efficiency rises to $\sim$ 35% by $10^{18}$ eV and peaks near $10^{20}$ eV at between 50-60%, depending on the specific station. Efficiencies for all stations and configurations are provided in additional Figures in App. A.2.

In this section, we describe the systematic uncertainties considered in the analysis. The impact of these systematics on $[A\Omega]_{\textrm{eff}}$ are shown in Fig. 9, and a table summarizing the magnitude of their effects at $10^{18}$ eV is provided in Tab. 1. We consider systematic uncertainties broadly in two classes. The first class is associated with theoretical uncertainties surrounding the neutrino-nucleon cross section and Askaryan emission, and are shown in Fig. 9 as solid bands, reported at the trigger level. The second class is associated with uncertainties in our understanding of the detection medium and our instrument. The latter are taken into account in setting the final limit as described in App. A, and are shown as dashed/dotted lines in Fig. 9 at the analysis level.

For the neutrino-nucleon cross section ($\sigma_{\nu-\textrm{N}}$), AraSim uses the model derived by Connolly, Thorne, and Waters (CTW) [18]. The upper and lower bounds for $\sigma_{\nu-N}$ are substituted for the central value in the simulation to estimate the effect of the uncertainty on the simulated $[A\Omega]_{\textrm{eff}}$ at the trigger level. In the CTW model, the uncertainties on $\sigma_{\nu-\textrm{N}}$ are large and grow as a function of energy, exceeding 100% above $10^{21}$ eV. At $10^{18}$ eV the uncertainties on the trigger-level effective area due to the cross-section are estimated at -15%/+18%. In Fig. 9, for comparison we also show the uncertainties if we use an alternative cross-section developed by Cooper-Sarkar et. al.(CS) [19] which has smaller uncertainties at high energies by about a factor of four. We additionally studied $d[A\Omega]_{\textrm{eff}}/d[\sigma_{\nu-\textrm{N}}]$, and find it to be approximately linear; for example, at 1 EeV, a 10% increase in $\sigma_{\nu-\textrm{N}}$ corresponded to a 10% increase in $[A\Omega]_{\textrm{eff}}$.

For the Askaryan emission, AraSim implements a modified version of the model derived by Alvarez-Muñiz et. al. [48]. A full description of modifications is provided elsewhere [31], but the primary differences arise due to the inclusion of of the LPM effect by Alvarez-Muñiz but not by AraSim, and in AraSim’s use of functional parameterizations for the shower profile instead of directly simulated shower profiles. The relative difference between waveform amplitudes produced by AraSim, and those derived from a full shower Monte-Carlo are at most ${\sim}12$% [49]. We conservatively estimate the effect of this systematic uncertainty by reducing or increasing the simulated field amplitude by $\pm 12\%$ and assessing the change in $[A\Omega]_{\textrm{eff}}$ at the trigger level. The relative difference between the default parameterization and the scaled parameterization has a maximum value of about 25% near $10^{16}$ eV, and starts falling as energy increases. This is because at high energies the instrument acceptance becomes dominated by geometric effects (ray tracing, etc.) and not signal amplitude. At $10^{18}$ eV the estimated uncertainties due to the Askaryan emission model are -11%/+13%.

In the second category of uncertainties, we consider those arising from our detector response and from measurements of quantities such as the index of refraction in ice and the attenuation length of radio waves in ice. These systematics are included in our calculation of the final limit. We consider uncertainties associated with (1) the attenuation length ($L_{\textrm{att}}$) of South Pole Ice and (2) the depth-dependent index of refraction ($n(z)$) of South Pole ice, (3) the calibration of the ARA signal chain, and (4) the triggering efficiency of the detector.

The model for the attenuation length ($L_{\textrm{att}}$) of South Pole ice was derived from data taken with the ARA Testbed prototype [26]. Confidence bands providing an upper and lower limit on $L_{\textrm{att}}$ are given in the model. To set upper/lower limits on our sensitivity, in AraSim, the upper and lower bounds for $L_{\textrm{att}}$ are substituted for the central value. At $10^{18}$ eV the uncertainty on the analysis level effective area due to uncertainties in attenuation length are -8%/+50%.

The model for the depth-dependent index of refraction $n(z)$ was obtained by fitting data obtained by the RICE experiment [41]. The data was fitted with an exponential as a function of (negative) depth $z$ of the form $n_{d}-(n_{d}-n_{s})e^{z\cdot n_{c}}$, finding the following parameter values and their respective uncertainties: $n_{d}=1.788\pm 0.016,$, $n_{s}=1.359\pm 0.022$ and $n_{c}=0.0132\pm 0.0017\,\textrm{m}^{-1}$. We recalculate the sensitivity, setting all parameters to their upper and lower limits simultaneously. The lower (upper) limit generally corresponds to a slower (faster) transition from surface to deep ice, and correspondingly have a smaller (larger) geometric acceptance for neutrinos. Additionally, since we do not change the ice-model assumption used to reconstruct the incoming direction of the RF emission as discussed in Sec. III.4, this systematic uncertainty also captures errors which may be present if the true ice model for radio wave propagation does not match that used for reconstruction. At $10^{18}$ eV the uncertainties on the analysis level effective area due to the index of refraction model are 5%.

We consider four sources of uncertainties that exist in the signal chain. They are the transmission coefficient $t$ representing the impedance mismatch between the ice and the antenna, as well as between the antenna and the coaxial cable, the ambient noise power received $N_{\textrm{ant}}$, the signal chain noise power $N_{\textrm{sc}}$, and the antenna directivity $D$. We follow the treatment used in the previous ARA result [32] where we consider the system signal-to-noise ratio representing the ratio of input signal power to total system noise power in a given channel:

with $P_{\textrm{sig}}$ being the received signal power. The four sources of uncertainty translate to an uncertainty in $SNR_{\textrm{sys}}$ by standard error propagation, which is then implemented as an uncertainty in the antenna gain $G$ in code ($\Delta G=\Delta SNR_{\textrm{sys}}/P_{\textrm{sig}}$). In line with previous ARA work, here we only consider the case where the effective gain of the instrument is reduced, providing a conservative estimate of our sensitivity. This is done because we lack sufficient calibration data at this time to constrain the upper bound on the gain. The VPol antenna gain has an overall estimated uncertainty of -10%, while the HPol antenna gain is estimated at -32%. The modified gain values are substituted in the simulation to assess the impact of this uncertainty, and the uncertainty at $10^{18}$ eV is found to be -3%.

For the systematic uncertainty associated with the trigger efficiency of the detector as a function of RPR, $\epsilon(RPR)$, we compare the simulated trigger efficiency $\epsilon_{\textrm{sim}}(RPR)$ to the measured trigger efficiency in calibration pulser data $\epsilon_{\textrm{dat}}(RPR)$: ${\Delta\epsilon=\epsilon_{\textrm{dat}}(RPR)-\epsilon_{\textrm{sim}}(RPR)}$. We measure $\epsilon_{\textrm{dat}}(RPR)$ by varying a tunable attenuator on the local calibration pulsers described in Sec. II and counting the number of calibration pulsers recorded. Using Arasim we find that the uncertainties on the trigger efficiency decreases the simulated $[A\Omega]_{\textrm{eff}}$ from between 2-5% depending on energy, and at $10^{18}$ eV the size of the effect is -3%.

We observed in previous calibration exercises that the stations trigger inefficiently on calibration pulsers whose direct ray-tracing solution intercepts the array at an angle steeper than $-25^{\circ}$ from from horizontal; this can be seen in Ref. [50], where there is a deficiency of triggers in A2 and A3 after the pulser is lowered below 1300m depth, despite the pulser being lowered to a total depth of 1700m. Therefore, for the calculation of $[A\Omega]_{\textrm{eff}}$ used in the limit, we conservatively exclude neutrino simulated events with the same ray-tracing conditions. This results in a ${\sim}10-30\%$ reduction in sensitivity, depending on energy. Excluding these steeply upgoing events is a conservative approach, as more exhaustive future studies might reveal that the cause of the trigger inefficiency to the calibration pulses does not have the same effect on neutrino events.

In this paper, we present constraints on the flux of UHE neutrinos between $10^{16}$ and $10^{21}$ eV from four years of data in A2 and A3. We have presented a description of the livetime and the instrument, and detailed the cuts used to eliminate backgrounds in two complementary, blind analyses. The resultant limit from this search is the strongest limit set by ARA to date, and the strongest limit set by an in-ice radio neutrino detector above $10^{17}$ eV. The result utilizes more than quadruple the livetime of the previously published ARA analysis, and maintains reasonable efficiency to neutrinos while remaining general to signal shape and not requiring costly cuts on livetime in Austral summer or angular cuts in the direction of anthropogenic sources like South Pole Station. We are encouraged that the two analyses, which leveraged complimentary sets of reconstruction and analysis tools, have similar sensitivity and produced consistent expected limits within 15% for all energy bins.

Post-unblinding, we were additionally able to further study our surface related backgrounds. As discussed previously, we observed zero events, consistent with our background estimates, including our estimated $10^{-3}$ events from above the surface. If we check the data taking runs that were excluded pre-unblinding because of the presence of large amounts of surface noise, we do observe a few events passing all cuts. We additionally are able to roughly estimate the probability of a single event of surface-origin being misreconstructed as coming from within the ice. To do so, we take the product of the fraction of runs in which we observe only one surface event and multiply by our estimated misreconstruction rate. We estimate the misreconstruction rate by taking the ratio of the number of events reconstructing inside the ice, relative to those reconstructing outside the ice, in the surface noisy runs. For example, in A3, we find that there may be approximately 0.2 such “misreconstructing singlets.” We note that this estimate is biased to larger values, because in order to measure the misreconstruction rate, we rely on the number of events reconstructing inside the ice in runs which demonstrate large amounts of surface noise, and were decided pre-unblinding to be unfit for analysis. These two post-unblinding studies demonstrate the role of the surface-noisy cut in the present analysis, and represent an opportunity for growth during the development of future reconstruction techniques.

We underscore several important features of this newest result. First, it demonstrates ARA’s capability to analyze its growing dataset. Compared to our previous result, which analyzed the first 10 months of data from stations A2 and A3 [32], this analysis leverages data from four years of data-taking in each of the two stations. After removing intermittent periods of downtime we have about 1100 days (75%) of livetime that was good for analysis for each station. This amounts to 2162 days of combined livetime. This analysis is therefore the first ARA result to analyze $\mathcal{O}(10)$ station-years of data. This demonstrates the capability to analyze our growing dataset, which will be important as ARA looks to the future. There is roughly 4080 additional days of livetime awaiting analysis on archive, with the analysis pending ongoing calibration efforts. With the full five-station ARA array collecting data since January 2018, the data set is expected to roughly double again by 2022 (total of approximately 11k days of livetime). In Fig. 8 we additionally show the projected trigger-level single-event sensitivity that the five-station ARA5 array can achieve with data that will have been accumulated through 2022. As can be seen, ARA is poised to be the leading UHE neutrino detector above $10^{17}$ eV; the IceCube and Auger experiments will also accumulate additional livetime amounting to about 40% and 25% increases over their respective published limits.

Second, the analysis maintains reasonable efficiency (${\sim}35$% at $10^{18}$ eV, and reaching 50% efficiency near a voltage signal-to-noise ratio (SNR) of 6) while remaining general and not relying on quantities that are strongly model-dependent, such as a correlation with a signal template. This is advantageous because although the Askaryan signal has been observed in the laboratory and in the atmosphere from cosmic-ray air showers [51, 52, 53], it has never before been observed in a dense media in nature.

In line with our previous two-station result [32], this analysis did not require excluding data recorded during the Austral summer, nor did it require geometric rejection regions specifically in the direction of the South Pole. In the prior analyses of the Prototype station [31, 38], 31% of livetime was lost due to anthropogenic activities during the Austral summer, as well as 9% due to the detector’s solid angle coverage in directions near the South Pole.

We note three challenges overcome in these analyses that have resulted in improvements moving forward, especially as the ARA dataset continues to grow, the diversity of the array increases, and the field looks forward to a large-scale radio array in IceCube-Gen2. The first challenge was managing the time-dependent nature of the ARA instruments. Some of the time dependent nature is owed to the different data taking configurations, as described in App. A.2, and often reflect improved understanding of the instrument and the ice. For example, an early trigger configuration led to triggering signals being off-center in the digitized waveform, and this was later corrected. The change to the readout length was made after Monte Carlo studies revealed that longer readout windows increased the probability with which a station records both the direct and refracted/reflected pulse that are possible because of the depth dependent index of refraction. Since learning from these processes, we have reached more stable operations configurations, and are working on additional streamlining. Some time dependence is owed to changing detector characteristics; for example, for some periods of time in ARA station 3, a digitization board exhibited a high amount of readout noise. Such time-dependent detector features required adjustments to analysis algorithms and analysis thresholds. As a result of the analyses described herein, identification of such time periods has also been considerably streamlined.

The second challenge was improvement in intra-collaboration communication between the ARA operations and analysis teams. In many cases, periods of livetime that were contaminated with calibration activity were recorded in operations reports, but were only later accounted for the analyses. We plan to work to streamline this pipeline for future ARA analyses.

The third challenge was managing anthropogenic activity from the South Pole over several Austral summers. Despite most human activity being isolated nearly two miles away, the analysis requires aggressive cuts on downgoing signals, which eliminated 10-30% of neutrino events. Improvements to reconstruction algorithms to more confidently reject downgoing events without requiring such substantial cuts on solid angle, or to more confidently reconstruct events with low hit-multiplicity, will improve the analysis efficiencies in the future.

The main authors of this manuscript were Brian Clark, Ming-Yuan Lu, and Jorge Torres, with Brian Clark and Ming-Yuan Lu leading the data analysis for this result. The ARA Collaboration designed, constructed and now operates the ARA detectors. Data processing and calibration, Monte Carlo simulations of the detector and of theoretical models and data analyses were performed by a large number of collaboration members, who also discussed and approved the scientific results presented here. We are grateful for contributions and discussions from Lucas Smith and Suren Gourapura.

We are thankful to the National Science Foundation (NSF) Office of Polar Programs and Physics Division for funding support through grants 1806923, 1404266, OPP-902483, OPP-1359535, and 1607555. We further thank the Taiwan National Science Councils Vanguard Program NSC 92-2628-M-002-09 and the Belgian F.R.S.-FNRS Grant 4.4508.01. We also thank the University of Wisconsin Alumni Research Foundation, the University of Maryland, and the Ohio State University for their support. B. A. Clark thanks the NSF for support through the Graduate Research Fellowship Program Award DGE-1343012 and the Astronomy and Astrophysics Postdoctoral Fellowship under Award 1903885, as well as the Institute for Cyber-Enabled Research at Michigan State University. A. Connolly thanks the NSF for CAREER Award 1255557 and Award GRT00049285 and also the Ohio Supercomputer Center. K. Hoffman likewise thanks the NSF for their support through CAREER award 0847658. S. A. Wissel thanks the NSF for support through CAREER Award 1752922 and the Bill and Linda Frost Fund at the California Polytechnic State University. A. Connolly, H. Landsman, and D. Besson thank the United States-Israel Binational Science Foundation for their support through Grant 2012077. A. Connolly, A. Karle, and J. Kelley thank the NSF for the support through BIGDATA Grant 1250720. D. Besson and A. Novikov acknowledge support from National Research Nuclear University MEPhi (Moscow Engineering Physics Institute). K. Hughes thanks the NSF for support through the Graduate Research Fellowship Program Award DGE-1746045. A. Vieregg thanks the Sloan Foundation and the Research Corporation for Science Advancement. R. Nichol thanks the Leverhulme Trust for their support. K.D. de Vries is supported by European Research Council under the EU-ropean Unions Horizon 2020 research and innovation program (grant agreement No 805486).

Finally, we are thankful to the Raytheon Polar Services Corporation, Lockheed Martin, and the Antarctic Support Contractor for field support and enabling our work on the harshest continent.