Reconstructing the arrival direction of neutrinos in deep in-ice radio detectors

Neutrinos with energies up to 10 PeV have been detected numerously and IceCube has discovered a component of astrophysical neutrinos above the atmospheric neutrino background [1]. These extraterrestrial neutrinos are expected to be created in extreme cosmic sources that accelerate charged particles, cosmic rays, to high energies, which produce secondary particles in their interactions with ambient matter and photon fields: neutrinos e.g. [2,3,4]. Several point sources are revealing themselves as source candidates for neutrinos, and therefore of the parent cosmic rays, either by an excess of neutrino events or through a coincident detection with multi-messenger particles [5,6,7]. However, the first neutrino source still needs to be identified at \(5 \sigma \) significance. As neutrinos typically carry \(\mathcal {O}(1/20)\) of the energy of the parent nucleon [8], a different yet undiscovered population of neutrinos will be the result of the interactions of ultra-high energy (UHE) (> EeV) cosmic rays, for which the energy spectrum is accurately established [9, 10]. Radio detectors are proposed for the detection of the resulting UHE neutrinos. Among these UHE neutrinos also cosmogenic neutrinos are expected, which result from the interaction of cosmic rays with photon fields during their propagation towards Earth [3].

The currently largest neutrino detectors use optical Cherenkov light stemming from the secondary particles that are created when the neutrino interacts (e.g. [1, 11]). Since this light has a propagation length of \(\mathcal {O}(100)\) m due to scattering and absorption in the relevant dense media [12], effective volumes are limited, as the volume is required to be densely instrumented. The secondary particle showers also emit nanosecond-duration radio pulses, which can travel for \(\mathcal {O}\)(1 km) in dielectric media such as ice. Therefore sparse radio arrays can be built to cover large volumes. Where in-ice radio detectors are in principle sensitive to the highest energies (>10 PeV), they rapidly lose sensitivity below 100 PeV, due to the irreducible thermal noise of the antennas and the surroundings, which prevents the detection of very low amplitude signals.

This work focuses on radio detectors using ice as the detection medium. The technique is used for in-ice radio neutrino detectors in the pilot arrays ARIANNA [13] and ARA [14], in RNO-G, currently being constructed in Greenland [15], and the proposed radio component for the next generation of IceCube, IceCube-Gen2 [16]. RNO-G is scheduled to be completed with 35 stations in 2026 and will be complementary to IceCube-Gen2 in the northern hemisphere. It will be the first in-ice radio array of sufficient scale to have a realistic chance of discovering the first radio neutrino. Such a discovery will be facilitated by a good angular resolution. Its impact for astrophysics will be also enhanced by a reconstruction of the neutrino arrival direction [17]. Besides aiding in the identification of neutrino sources, a good resolution in zenith helps to firmly establish the diffuse UHE neutrino flux, as not only the flux but also the neutrino cross section is not yet known at these energies [18]. IceCube-Gen2, with a planned effective volume of more than 1600 \(\hbox {km}^3\) sr, will reach the discovery space of almost all predictions of UHE neutrino fluxes, e.g. [19,20,21,22,23].

In-ice radio neutrino detectors typically consist of a cluster of in-ice antennas forming a self-triggering station, with distances between stations of \(\mathcal {O}(1)\) km such that there is little overlap in effective volume for individual stations. It is estimated that only 10% of the triggering neutrinos will be detectable in two stations. Radio neutrino stations are typically classified as shallow stations, deep stations, and hybrid stations. Shallow stations contain broadband antennas buried 1 m to 3 m below the ice surface. They are limited in effective volume, because the gradually changing refractive index with depth (due to the increasing ice density) results in parts of the ice from which the radio signal is not able to reach the detector. Deep stations contain deep in-ice holes (down to depths of 100–200 m) in which strings of antennas are installed. They reach down to the bulk ice, providing a large effective volume, but these holes are narrow, \(\mathcal {O}\)(30 cm), and therefore limit the geometry of the antennas. Hybrid stations are a combination of a deep and shallow station. Both RNO-G and IceCube-Gen2 will have hybrid stations.

1.1 Previous work

A number of reconstruction algorithms for different aspects of radio particle detectors have been developed and resolutions are highly dependent on the specific station design and detection medium.

The reconstruction of the signal arrival direction is best studied. It depends primarily on how well the antenna position and the instrument timing is known, as well as the assumed waveform model. In air showers, direction resolutions of better than 0.5\(^\circ \) are common, if many antennas detect the signal and a hyperbolic waveform fit can be used, as shown for e.g. LOPES [24] and LOFAR [25]. For in-ice experiments, the ARIANNA experiment (shallow antenna installation) reported an angular resolution of \(<1^\circ \) for a pulser lowered into the ice [26]. The ARA collaboration (deep antennas) has reported a similar angular resolution for pulsers installed at large depths [14].

Some aspects of reconstructing the signal polarization are also relatively well studied. Air shower experiments have shown that the statistical uncertainty on the angle of polarization is inversely proportional to the signal-to-noise ratio approaching \(0.1^{\circ }\) for clear signals [27]. Additional systematic uncertainties may stem from the antenna modeling, but were estimated to reach better than 10% [28]. The reconstructions even allowed for the measurement of the small circular component in the in-air emission, stemming from the time difference between geomagnetic and Askaryan emission [29].

The situation in ice is more complex for two reasons: Glacial ice shows a gradient in the index of refraction that significantly curves trajectories and the ice is known to exhibit birefringent properties that will influence how well the polarization can be reconstructed [30, 31]. So far, reconstructions have taken into account the bending, but mostly ignored the birefringence. Second, in the reconstruction of the electric field, antennas at different positions have to be combined rather than using dual-polarized antennas like in air, since antenna sizes are restricted due to installation constraints. This adds additional timing and positioning uncertainty.

For shallow antennas a polarization resolution of order 1\(^\circ \) was obtained based on cosmic-ray data, propagating only briefly through the ice, but in excellent agreement with simulations [32]. ARIANNA also reported a resolution of 3\(^\circ \) for in-ice pulser data, including emitter uncertainties, albeit only probing a small range of polarization angles with high signal-to-noise data [26].

The ARA experiment has reported on a number of indications for birefringence [33,34,35]. Their data probe much longer distances and more polarization angles, but a complete model of the signal propagation of neutrino pulses is still outstanding, preventing the full electric-field reconstruction for the polarization reconstruction of pulser signals.

Aside from experimental evidence, several simulation studies have been performed to study reconstruction capabilities of in-ice radio stations. Machine learning algorithms show a lot of promise for neutrino angular reconstruction. For shallow antennas, a space angle difference of roughly \(4^{\circ }\) has been found [36]. Similar results have been reported for ARA, albeit less studied in detail at this point [37].

The neutrino energy resolution for deep antennas has been studied in [38]. The presented method uses Information Field Theory (IFT) [39] to reconstruct the electric field at the antennas. The amplitude of the measured electric field pulses scale linearly with the energy and inversely with the distance to the interaction vertex. Therefore, the energy reconstruction is highly dependent on the distance to the interaction point of the neutrino, i.e. the vertex position.

Due to the in general low amplitude of the signal pulses, standard methods that unfold the voltage data with the detector response to obtain the electric field typically overestimate certain signal parameters and perform worse than methods where noiseless waveforms are matched with the data. Examples are the IFT approach and the forward-folding approach, where noiseless voltage waveform expectations are obtained by forward-folding an electric field assumption through the ice and detector response [40]. Forward-folding will also be used as the basis of this work. It has previously been applied to shallow antennas, resulting in an angular resolution of 3\(^\circ \) [41].

Besides the applicability to low amplitude events, the model-dependency of our method has the advantage (as compared to the current implementation of IFT) that it takes into account the different electric fields at different antenna positions, i.e. it knows where the antennas are.

After a didactic overview of neutrino angular reconstruction in Sect. 2, explaining the complexity of pinpointing the neutrino arrival direction due to the shower geometry, we will describe the algorithm and reconstruction approach in Sect. 3. We show the performance of the algorithm and the resulting angular resolution in Sect. 4. We will also discuss the impact of different selection cuts on the obtained resolutions. Furthermore, we will show event contours of typical events, as well as the point spread function, and discuss the implications for point-source studies. Our conclusions will be summarized in Sect. 5.

In-ice radio detectors measure the radio emission stemming from particle showers induced by neutrinos interacting in vast ice volumes. The radio emission comes from an accumulated negative charge excess which varies over time due to the shower development, the Askaryan effect [42]. It results in coherent nanosecond long radio pulses, emitted on a cone with its maximum amplitude at \(\approx \) 56\(^\circ \) (for deep ice in Greenland), the Cherenkov angle \(\varTheta _{C}\).

Illustration of the shower geometry and the radio cone on which the signal is typically detectable (right, not to scale), as well as the impact on the reconstruction of a single event (top left) and the point spread function (bottom left). The point spread function is shown based on events with a number of different viewing angle and polarization combinations to illustrate the difference with the single event reconstruction. The radio signal is observed at the detector under the viewing angle, i.e. the angle between the shower axis and the radio signal, which is typically a few degrees off from the Cherenkov angle. The electric field is polarized towards the shower axis due to the Askaryan effect. A measurement of the signal direction, viewing angle and polarization angle results in the reconstruction of a single neutrino as illustrated in the top left figure. Multiple neutrinos from a point source are reconstructed as shown in the bottom left figure

Full size image

Due to the cone, the signal emission direction is not equal to the neutrino arrival direction. In addition, due to the varying index of refraction, the direction of the signal changes during propagation such that the emission direction is not the same as the direction when the signal arrives at the detector. However, a measure of the signal direction, the angle between the radio emission and the neutrino direction (the viewing angle), and the polarization of the electric field are sufficient to fully constrain the neutrino direction. This leads to a complex number of steps involved to reconstruct the neutrino direction, which will be elaborated in the following sections.

2.1 Connection of neutrino direction and signal parameters

As the particle shower propagates at a speed close to c, and hence faster than the speed \(v = \frac{c}{n}\) of the electromagnetic waves, with n denoting the refractive index, radio waves arrive simultaneously at an observer at \(\varTheta _{C}\) from the shower axis. This observer sees fully coherent radio waves. For slightly off-cone angles, the coherence for high frequencies is lost, resulting in a lower cut-off frequency and rapidly decreasing signal amplitude. Therefore, the electromagnetic waves can only be observed on a cone around the shower axis with opening angles close to the Cherenkov angle \(\varTheta _{C}\pm 10^\circ \) for detectors sensitive to MHz frequencies.

It also means that the frequency content of the observed pulses contains information on the observed viewing angle of the radio emission with respect to the shower axis. The negative charge excess in the shower front, leads to a polarized electric field with its electric field vector pointing towards the shower axis. Therefore, combined measurements of the signal emission direction, viewing angle, and polarization suffice to reconstruct the direction of the original neutrinos. This is illustrated in Fig. 1, left. In reality, the method reconstructs the direction of the particle showers. However, the shower is very strongly forward-boosted already at 1 TeV [43], such that at the energies relevant here the angular deviation between the shower and the original neutrino is completely negligible.

2.2 Modeling of neutrino radio signals

Radio emission from showers can be calculated from first principles in microscopic particle simulations, e.g. [44, 45]. This is, however, computationally expensive and only needed, if signal parameter accuracies below 10% are needed. For neutrino detection, where large uncertainties stem from the propagation in ice, parameterizations are currently sufficient. Modeling the electric field expected from the Askaryan emission from a particle shower is best done using a semi-analytical model, that calculates the time domain waveform directly from the charge-excess distribution in the shower via convolution with a form factor that only depends on the shower type, i.e. hadronic and electromagnetic showers [46]. The model shows good agreement (\(3\%\)) with Monte Carlo codes. At the cost of a loss of accuracy, simpler (and thus faster) parameterizations are also available.

In addition to the parameters discussed in Sect. 2.1, the Askaryan emission also depends on the type of particle shower initiated by the neutrino. For \(\nu _\mu \), \(\nu _\tau \), and neutral current interactions from \(\nu _e\), only a hadronic shower is created, but for \(\nu _e\) charged current interactions, the electron additionally gives rise to an electromagnetic shower. At ultra-high energies, electrons are subject to the LPM effect, resulting in a more irregularly shaped shower profile than hadronic showers [47, 48], and the presence of multiple showers can cause interference at the detector. Accordingly, the reconstruction of these events is expected to be more challenging. The algorithm outlined in this article was therefore designed for hadronic events. We do, however, include all events in our simulation, and discuss the impact on performance in Sect. 4.

The current state-of-the-art neutrino simulation code, NuRadioMC [49], provides an interface to most of the models. NuRadioMC modularizes the simulation problem into event generation (neutrino vertices), signal generation (chosen from a parameterization), signal propagation, and signal detection (flexibly including various detector types). NuRadioMC is open source and under continuous development, adding for example more complex propagation models and better descriptions of the detectors. All simulations described here are performed with NuRadioMC.

2.3 Signal propagation

NuRadioMC uses the decoupling of signal propagation and signal generation, i.e. for a given neutrino interaction position the path traveled through the detector is calculated using ray tracing. The electric field as propagated on this path and observed at the antenna is then calculated. The ray tracing approximation holds when the wavelength is much larger than the relevant features in the ice. In most cases this is sufficient, i.e. the refractive index is slowly increasing and the ice is homogeneous.

The refractive index of polar ice is increasing with depth, reaching a constant value of 1.78 for deep ice at \(\approx \) 100 m (Greenland). The resulting changing speed of the radio waves causes the waves to bend towards regions of higher refractive index, so that signals from neutrino vertices below an observer can also arrive refracted from above. At an air-ice, ice-water or ice-rock interface, radio waves are (partly) reflected. During propagation, the radio waves are attenuated, the attenuation length of the medium is both depth and frequency dependent. All of these features are taken into account in NuRadioMC, resulting in so-called direct, refracted, and reflected signal propagation paths. Horizontal propagation of signals near the ice-air surface or due to density fluctuations in the firn [50, 51] are currently not taken into account.

Recent data has shown that polar ice can be birefringent, i.e. the ice crystal orientation causes the refractive index to be dependent on the direction and polarization of the radio waves which results in a time-delay for different polarization components. Work is ongoing in understanding and simulating this effect e.g. [30, 31, 52], but the modeling is currently not mature enough to include this in our simulation study. Furthermore, layers in the ice with density fluctuations can cause partial reflection of the radio waves. These layers are observed at South Pole [35] and Greenland [53]. Also their effect on the angular reconstructing has to be studied in future work.

2.4 Detector response

Deep in-ice stations typically contain two kinds of antennas to measure orthogonal electric field components, designed to be predominantly sensitive to either the vertical polarization (Vpol) or the horizontal polarization (Hpol). The Vpol is designed to be high gain and broadband in the \(\mathcal {O}\)(100–600) MHz range, such that information on the viewing angle can be extracted from the frequency content. The Hpol design, typically lower in gain due to the geometry constraints by the borehole, is designed to have sufficient overlap in the frequency content with the Vpol to aid polarization reconstruction.

In-ice neutrino radio arrays are designed to optimize effective volumes, hence the independent triggering stations are located on a sparse grid. The bulk of neutrinos will trigger only a single station, therefore each station is designed to be capable of reconstructing neutrino properties. A schematic of the deep station used in this work is shown in Fig. 1 (right). It contains three strings of antennas at lateral distances of 35 m. The holes are drilled to such a depth that the deepest point of the strings is below the firn, to maximize the effective volume. At the bottom of the power string closely spaced Vpols (1 m) are located. By beamforming the waveforms in these antennas, which is an effective noise reduction, they function as a low-threshold trigger called the phased array. Furthermore, the power string is equipped with Hpols located right above the phased array, and Vpols at larger distances. The other two helper strings have Vpols and an Hpol located at the lower end of the string.

Simulations take into account both the antenna locations, as well as their complex performance (gain and group delay) to incoming signals. Also, the effect of the electronics chain (amplifiers, cable losses, digitization, triggering) is simulated, to be able to study the reconstruction efficiency. The station modeled here closely resembles an RNO-G station as currently under construction.

2.5 Event contour and point spread function

Due to the number of steps involved in reconstructing the neutrino arrival direction for radio neutrino detectors, there is a distinct difference between the shape of a single event contour and the point spread function (PSF). As will be elaborated further in detail, both should only be approximated in two dimensions. This is illustrated in Fig. 1.

A measure of the signal arrival direction restricts the neutrino direction to a broad band, because the amplitude at the detector will only be large enough to be detected for angles a few degrees off the Cherenkov cone. A viewing angle measurement narrows this band. An additional polarization measurement restricts the location of this band. This results in ellipse like single-event uncertainties for good quality events.

Reconstructing the polarization is most challenging for deep in-ice radio detectors. The design of the antennas is constrained by the borehole diameter, typically \(\approx \) 30 cm wide, due to the available drills. This makes designing antennas sensitive to the horizontal signal component (Hpol antennas) more challenging than the ones sensitive to the vertical component (Vpol antennas). This then results in a significantly smaller gain for the Hpols, with them often measuring no detectable signal. Therefore, the main contribution to the angular resolution is the uncertainty in the polarization measurement resulting in ellipse-like uncertainty contours with very large axis ratios.

A neutrino source observed with an in-ice radio station will be seen smeared out due to the reconstruction uncertainties, which results in the sum of ellipse-shape contours. Due to the geometry of the shower and the corresponding radio-cone, and the fact that the Vpol antennas are used for triggering, the range of interaction positions of the neutrino in the ice (thus signal arrival directions) that can trigger a station is limited. Since the orientation of the ellipse is determined by the signal arrival direction, this results in orientations of the ellipse that are forbidden, resulting in a ’bow-tie’-shaped PSF.

It has been first shown in [40] that forward-folding is a suitable reconstruction tool for radio detectors. In the forward-folding approach, an analytic description of the electric field is forward-folded through the propagation and system response to obtain waveforms of the voltage received at every antenna. These waveforms are compared with the recorded voltage data and a test statistic is minimized to obtain the best fitting signal parameters. This approach has shown to work better for reconstructing low signal amplitude radio pulses than unfolding the detector response. Due to the contribution of noise in the voltage data, unfolding the (noiseless) detector response systematically overestimates the signal contribution when the detector response expects a low value.

For reconstructing the neutrino direction, the forward-folding approach of [40] has been extended to use antennas with large distances to each other, where it no longer holds that the same signal is measured. This extension has the benefit that using more antennas in the same reconstruction effectively reduces the noise contribution and therefore improves the obtainable resolution. Furthermore, combining antennas across the entire station results in a more accurate measurement of the polarization due to the vertical and horizontal spatial separation than obtainable by only using nearby antennas.

3.1 Electric field model

For reconstructing the radio pulses an analytical description of the electric field as provided by [54] is used which gives a parameterization of the electric field in the frequency domain, where the amplitude is linearly rising with frequency up to a cut-off frequency (\(\approx \) 1.2 GHz for on cone angles), which reduces for off-cone angles. The parameterization depends on the viewing angle assuming the emission is coming from the shower maximum, i.e. the point which most of the emission originates from. Besides the viewing angle, the electric field parameterization depends only on the shower energy \(E_{\text {sh}}\) which determines the signal amplitude. Since we use a parameterization in the frequency domain, signals coming from different parts of the shower arrive at different times, i.e. the phase, is not taken into account. Consequences for the reconstruction procedure are explained in Sect. 3.3. It should be noted that the parameterization used to reconstruct the data, is simpler than the one used to model that data (see Sect. 2.2), which should make the approach more robust towards slight mismodeling of the signal.

Since the electric field model depends on the viewing angle, the interaction point of the neutrino is required to compute the viewing angle under which the radio emission is seen at the antenna for a given neutrino direction. Therefore, reconstructing the interaction point is the first step in the reconstruction procedure.

3.2 Neutrino vertex reconstruction

We use the vertex reconstruction algorithm, which effectively reconstructs the shower maximum, as described in [38]. As there have been no major changes to the algorithm, we only provide a brief summary below, and refer to [38] for a more detailed description.

The relative pulse arrival times in the different antennas depend on the vertex position. Therefore, the vertex position can be triangulated by shifting the voltage traces depending on the hypothesized vertex position, and maximizing the resulting antenna-to-antenna pairwise correlation. In fact, rather than correlating the voltage traces directly with each other, they are first correlated with a neutrino template, reducing the effect of accidental correlations in the noise in each antenna. In addition, a bandpass filter with a passband of 95 MHz to 300 MHz is used before performing the correlation. The pulse shape in this band is determined mostly by the antenna and amplifier response, and varies only slightly depending on the arrival direction and viewing angle.

One additional challenge with respect to traditional triangulation arises due to the presence of multiple ray paths between vertex and antenna position, as explained in Sect. 2.3. This is taken into account by testing the time shifts for each possible ray type (’direct’ or ’refracted’/’reflected’) for each antenna. Finally, one additional subtlety arises due to the fact that the emission does not originate from the neutrino interaction vertex, but from the shower, whose core extends over the order of tens of meters (or more, for electromagnetic showers affected by the LPM effect [47, 48]). Therefore effectively the position of the shower maximum is reconstructed and not the vertex, which has been accounted for in the following steps.

3.3 Reconstruction algorithm

Although the reconstructed vertex position provides the expected relative time delays between the different antennas, the absolute pulse arrival time is needed in addition in order to determine the pulse windows to include in the fit. Pulse windows are used to limit the options for accidental correlations. The absolute time is determined by the ray-type selection algorithm. As in the case of the vertex reconstruction algorithm, a template is correlated with the voltage traces for each ray-type hypothesis, and the dominant ray type and its arrival time are determined by looking at the maximum total resulting correlation. This is illustrated in Fig. 2.

The pulse selection procedure. Left, top: the total template correlation for the phased array antennas for each ray-type hypothesis (’direct’ or ’refracted’/’reflected’), indicating here a better correlation for refracted. Left, bottom: The individual as well as coherently-summed (beamformed) pulses for each ray type. Right: Window selection for both ray types. The pulse windows to include in the fit are selected based on their expected locations relative to the pulse arrival times in the phased array, and a cut on the peak-to-peak SNR. Included pulse windows are highlighted in green, excluded ones in gray. The pulses corresponding to the ’phased array cluster’ (shown in boldface), consisting of the phased array and adjacent Hpol antennas, are always included in the fit

Full size image

Unlike in the vertex reconstruction algorithm, here only the four deep Vpol antennas in the phased array are used. As they form the trigger of the experiment, the likelihood of observing a clear (when beamformed) pulse here is highest. In addition, their proximity to each other reduces the relative timing error due to differences in viewing angle or uncertainties in e.g. the ice model.

Once the pulse arrival time and the ray type in the phased array have been identified, the approximate pulse windows for the other antennas for each ray type can be calculated. The exact pulse arrival time in the other antennas is however not known: first of all, the electric field model used in the reconstruction does not contain phase information, leading to a small timing shift in the antennas at larger baselines due to relative shifts in the viewing angle. Additional timing shifts are caused by the error in the vertex reconstruction or uncertainties in the antenna positions and ice model.

Therefore, in order to identify the exact arrival time within each pulse window, another correlation between the reconstructed waveform and the data is performed. For low signal amplitudes, this results in a position determined by thermal noise fluctuations, which results in an overestimation of the power in these pulses. To avoid this, only pulses which exceed a signal-to-noise (SNR) ratio of at least 3.5, defined as half the maximum peak-to-peak voltage (within the pulse window) divided by the root-mean-squared noise voltage \(\sigma _\textrm{noise}\), are included in the fit. This is illustrated in Fig. 2, right. Note that the pulses in the four deep Vpols of the phased array are always included, as well as the pulses of the two Hpols (H7,H8) directly above them, as these are close enough to allow the pulse arrival times to be determined directly from the phased array. Finally, antennas where the two expected ray types have overlapping pulse windows, causing the pulses to interfere, are excluded from the fit.

The \(\chi ^2\) statistic to be minimized is obtained by forward-folding the model prediction of noiseless signal waveforms with the appropriate in-ice propagation and detector response effects, and comparing this with signal data (or simulations), which include a contribution of thermal noise. If the noise is assumed to be approximately Gaussian, the \(\chi ^2\)-statistic is given by

with \(x_{i}\) the voltage data at sample i, \(f_{i}\) the model prediction at sample i, \(\theta _\textrm{view}\) the viewing angle, \(\phi _\textrm{pol}\) the polarization angle and \(E_\textrm{sh}\) the shower energy. Only the \(n_{\text {pulses}}\) that pass the SNR > 3.5 criteria (with a minimum of \(n_{\text {pulses}}\) = 6), with a window of \(n_{\text {samples}}\) are included.

The pulses in the Hpol are typically shorter than the pulses in the Vpol, and therefore a smaller fitting time-window \(\varDelta t\) is used (30 ns) than for the Vpol (60 ns). This results in a different number of \(n_{\text {samples}}\) (sampling rate \(\cdot \varDelta t\)) depending on the antenna type. \(\sigma _\textrm{noise}\) is the noise root-mean-squared of the antennas after filters, which are applied to reduce the noise contribution in regions where the antenna is not sensitive. A low-pass filter of 700 MHz and a high-pass filter of 50 MHz are used, resulting in a typical \(\sigma _\textrm{noise}\) = 11.6 mV for this particular set-up.

The reconstructed pulses for one Vpol (top) and one Hpol (bottom) antenna in an example event, compared to the noiseless (sim) signal. The non-shaded region in the waveform (left) indicates which part of the waveform is included in the fit. The shape of the frequency spectrum is dominated by the system response

Full size image

The results of the fitting procedure for one example event and two of the included pulses are shown in Fig. 3 and compared to the true (simulated) neutrino signal. Note that the result here is additionally constrained by the fit to the 10 other pulses (not shown).

In order to test the reconstruction algorithm, we generate a representative set of neutrino signals as they would be detected in a deep radio station. We reconstruct the entire set and discuss the reconstruction performance for individual events, several high quality subsets, and the implications for radio detectors.

4.1 Simulation set-up

We use NuRadioMC to generate \(\sim 6.5 \cdot 10^7\) neutrinos with energies ranging from \(10^{16} - 10^{22}\) eV, following a spectrum that combines an IceCube-like flux [55] with an additional contribution from cosmogenic neutrinos as defined in [56]. The neutrino arrival directions are generated randomly, although their attenuation due to propagation in the Earth is taken into account by weighting the resulting events appropriately.

The cumulative distribution of the obtained resolution for different aspects of the reconstruction and different subsets of events. Note the different scales on the x-axis. Shaded regions indicate the approximate \(1 \sigma \) statistical uncertainty

Full size image

The simulated detector is based on the design of the RNO-G station, as described in [15]. The power string and the helper strings (see Sect. 2.4) are separated horizontally by 35 m. The power string consists of a phased array containing four Vpol antennas at \(\sim \) 100 m depth, with two Hpol antennas directly above, as well as 3 additional Vpol antennas at depths of 80, 60 and 40 m. The two helper strings feature two Vpol antennas with one additional Hpol antenna directly above, at a depth of 100 m. The antenna responses for the Hpol and Vpol are used as available in NuRadioMC. The biconal Vpol is sensitive between 100 and 600 MHz with a resonant frequency of about 150 MHz. The Hpol is significantly lower in gain and has its highest response between 200 and 500 MHz. Both have a symmetric response in azimuth and are most sensitive at a zenith of 90\(^\circ \). Both are connected to an amplifier (Low-Noise Amplifier, [15, Figure 17]) that boosts the signal strength (almost) frequency independent up to 800 MHz.

We model the phased array trigger (four Vpols, 1 m spacing) with a \(2 \sigma \) threshold trigger for a single antenna, which is studied to be the equivalent for the trigger used by RNO-G as shown in [15, Figure 18]. A sampling rate of 2.4 GHz is used, similar to the RNO-G design sampling rate. We simulate infinite trace lengths (trace lengths are adjusted such that all arriving pulses are stored). The influence of shorter recorded waveforms is not studied in this work. While no significant changes are expected for, e.g. the RNO-G hardware (record length of 850 ns), the impact could be larger for deeper stations as planned for IceCube-Gen2 Radio [16].

The noise due to thermal fluctuations of the electrons in the antennas can be described with

with \(k_b\) the Boltzmann constant, \(\varDelta f\) the frequency band, and \(T_{\text {noise}}\) the noise temperature. We simulate the noise with \(T_{\text {noise}}\)= 300 K, corresponding to a \(\sigma _\textrm{noise} = 11.6\ \textrm{mV}\) root-mean-squared of the noise within our frequency band of 50–700 MHz. We simulate the noise in the frequency domain with amplitudes that are Rayleigh distributed, resulting in (almost) Gaussian distributed noise in the time-domain.

The electric field model of [46] for the Askaryan emission is used, as described in Sect. 2.2. The ice refractive index is modeled with a depth-dependent refractive index model, described with an exponential function, as suggested by [50] and a depth and frequency dependent attenuation model from [57].

Note that the simulations performed are a statistical realization in two ways, i.e. 1) the thermal noise fluctuations and 2) the electric field generated from statistically varying shower profiles.

With this neutrino flux and detector configuration, we obtain an integrated weight of 4273 triggering events, of which 1881 contain only a hadronic shower. The reconstruction algorithm is run for each of these, and the results are described in the following section.

4.2 Reconstruction performance

The results of the reconstruction are shown in Fig. 4. We distinguish the signal direction resolution, i.e. the angle between the true and reconstructed emitted signals, the viewing angle, polarization angle, and space angle, i.e. the angle between true and reconstructed neutrino directions. Of these, the first depends only on the performance of the reconstruction of the shower maximum. We show results both for hadronic events (solid lines), for which our code was originally developed, as well as those for all events (including electromagnetic interactions, dashed lines). Performance for this latter category is expected to be worse, due to the more irregular shape of electromagnetic showers experiencing the LPM effect, as well as potential interference between the emission from the electromagnetic and hadronic showers. Particularly for these events, it is expected that the performance could be improved by developing a dedicated algorithm and/or classifier, as has for example been investigated in a machine-learning approach in [36].

In addition to the results for the full event set, we show three different subsets of events that emphasize the different features of the reconstruction. These are selected by imposing a minimum cut on the SNR in certain antennas, here defined as the maximum amplitude of the true (noiseless) signal divided by the true noise root-mean-squared amplitude \(\sigma _\textrm{noise}\). The three subsets are:

• The analysis cut: this includes events that have a signal with an SNR of at least 2.5 in one phased array antenna, as well as at least one other antenna in the power string. In addition, a signal with an SNR of at least 2 is required in each of the helper strings. This cut primarily improves the quality of the shower-maximum reconstruction, which requires a signal in at least these four antennas (groups) to uniquely determine the origin of the emission.

• The quality cut (V): this includes events that pass the analysis cut, and in addition have a signal with an SNR of at least 3.5 in one of the two uppermost Vpol antennas. This cut aims to further improve the quality of the shower-maximum reconstruction, by ensuring the presence of a clear signal at at least 5 different baselines.

• The quality cut (H): this includes events that pass the analysis cut, and in addition have an SNR of at least 3 in at least one of the Hpol antennas. This leads to a significantly stronger constraint on the polarization angle, which is the dominant uncertainty in the space angle.

The percentage of events retained by each of these cuts is shown in Table 1.

The obtained resolution is limited by the shower-maximum reconstruction. When this step is successful, the signal direction is reconstructed with a median resolution of 0.1–0.2\(^\circ \). However, in particular for low-SNR signals, or shallow vertices, which experience more significant ray bending in the firn, there are long tails in the distribution corresponding to misreconstructed vertices. This leads to wrong signal windows (i.e. the fit may not even include the actual pulses) and an incorrect geometry, which in general cause the rest of the reconstruction to fail, too.

If the shower maximum is reconstructed successfully, the dominant uncertainty is the polarization angle. The viewing angle is a priori constrained to lie within a couple of degrees of the Cherenkov angle, and can be reconstructed relatively well depending on the amplitude of the signal in the Vpol antennas. On the other hand, partly because of the relatively weaker response of the Hpol antennas, as described in Sect. 2.4, the resolution of the polarization angle is about an order of magnitude worse than that of the viewing angle. This leads to a space angle resolution that closely resembles the polarization angle resolution. As expected, the biggest improvement here is observed by applying a cut on the minimum SNR in the Hpol antennas (Quality cut (H)).

The difference between purely hadronic events and those including electromagnetic showers is most clearly visible in the signal direction and viewing angle resolutions. This is mostly due to two factors; firstly, the larger spread in the shapes of the individual showers for electromagnetic showers, and secondly, the current reconstruction algorithm uses an analytic model for the Askaryan spectrum of hadronic showers, which is slightly different from that for electromagnetic showers. A dedicated algorithm that additionally targets electromagnetic showers will likely improve on these results.

The space angle distribution for events passing the analysis cut, for different shower types and energies. For hadronic events, high-energy showers in general give higher amplitude signals, leading to an improving resolution as a function of energy

Full size image

The space angle distribution for events passing the analysis cut, for different shower types and zenith bins. More vertical showers correspond to shallower interaction vertices, due to the emission angle of approximately \(56^\circ \), which in turn are harder to reconstruct

Full size image

The resolution additionally varies with the deposited shower energy and the neutrino zenith. The energy dependence for the space angle for events passing the analysis cut is shown in Fig. 5. In general, even after the cuts are applied, higher energy showers are usually higher in amplitude in both Hpol and Vpol antennas, leading to an improving resolution as a function of shower energy. However, for electromagnetic showers, this effect is almost completely offset by the increased spread in the shower shape due to the LPM effect at higher energies. Furthermore, for lower energies a slightly better resolution in space angle for electromagnetic showers is obtained than for hadronic showers.

The dependence on the neutrino zenith angle is shown in Fig. 6, again for events passing the analysis cut only. In this case, the resolution improves as the neutrino becomes more inclined (larger zenith). There are two reasons for this. Firstly, more vertical showers (\(< 60^\circ \)) are geometrically constrained, by the Cherenkov angle of \(\sim 56 ^\circ \) and the fact that rays bend downwards in the firn, to lie closer to the surface, which leads to a more challenging shower-maximum reconstruction, as mentioned previously. In addition, as the polarization vector points towards the shower axis, more inclined neutrinos include a stronger horizontally polarized component, which improves the reconstruction of the polarization angle. These effects are similar for both hadronic and electromagnetic events.

Distributions of the reconstructed neutrino directions in two dimensions, rotated and projected into a coordinate frame where the polarization angle points along the x-axis and the viewing angle along the y-axis. Shown are individual events, as well as contour lines for different cuts, for all hadronic events

Full size image

Finally, it is important to keep in mind that the resolution, in terms of space angle, is almost completely dominated by the polarization angle, and therefore can not be straightforwardly interpreted as the relevant characteristic for e.g. point source searches. It is however possible to consider the reconstructed events in a rotated and projected reference frame, such that the (large) uncertainty in the direction of the polarization angle lies along the horizontal axis, and the (much smaller) uncertainty along the viewing angle direction along the vertical axis (c.f. Fig. 1, top left). The resulting asymmetric distribution is shown in Fig. 7. For the smaller quantiles, the contours are strongly elongated along the polarization direction, with the size of the contour decreasing with analysis and quality cuts. The larger quantiles, however, include events with poor shower-maximum reconstruction, for which the reconstruction error is more random, leading to more symmetric contours (e.g. the 90% contour for quality cut (H)).

We can then compute a resolution in square degrees, corresponding to the surface of an ellipse covering N% of events. This is done in Fig. 8. The secondary x-axis shows the corresponding 1D angle that would cover the same area for a symmetric contour (i.e. a circle). The median 2D resolution for hadronic events passing the analysis cut, for example, is approximately \(17\ \textrm{deg}^2\), corresponding to a 1D angle of \(2.4 ^\circ \), compared to a median space angle of \(4.9^\circ \). This figure also emphasizes the importance of improving the viewing angle reconstruction for electromagnetic events; although the resolution in terms of space angle is very similar for most events, including electromagnetic events increases the median spread in 2D by a factor of 2–3 in area.

The cumulative distribution of the obtained resolution, in terms of area in square degrees. The secondary axis indicates the angle that would correspond to the same surface angle in the case of a symmetric angular distribution. Shaded regions indicate the approximate \(1 \sigma \) statistical uncertainty

Full size image

4.3 Event contour and point spread function

Next, we study uncertainty contours for individual events. Figure 9 shows the \(68\%\) uncertainty region for an event that passed the analysis cuts. The contour is obtained as follows: on a grid in zenith and azimuth around the best fit point the test statistic is computed and a 68% contour of constant \(\varDelta \chi ^2 = \chi ^2 - \chi ^2_\text {min}\) is drawn. The appropriate value of \(\varDelta \chi ^2\) is calculated by reconstructing the same event multiple times with different realizations of the electric field and the thermal noise, and determining the 68th percentile of \(\chi ^2_{\text {true}} - \chi ^2_{\text {rec}}\), where \(\chi ^2_{\text {true}}\) is the test statistic corresponding to the true direction and best fitting shower energy. Thus, the contour can be interpreted such that 68% of the repeated experiments will be reconstructed within this contour. This is shown in the bottom plot of Fig. 9.

Top: The 68\(\%\)-contour for an event of average quality that passes the analysis cuts, resulting in a space angle of 4\(^\circ \). The uncertainty region is an elongated and narrow asymmetric area. Bottom: Distribution of the reconstructed direction for the same event under the influence of statistical uncertainties, e.g. shower profile and thermal noise fluctuations

Full size image

Same as Fig. 9 (top panel), but for a high-quality event with large contribution in the Hpol antennas

Full size image

An example contour for an event with a large amplitude in the Hpol is shown in Fig. 10. A \(68\%\) area of 4.5 \(\hbox {deg}^2\) is obtained. As can be seen from the figure, for small polarization angles the uncertainty region in zenith and azimuth can be closely approximated by an ellipse.

The event contours for neutrinos detected with the deep station, even those of high quality, are highly asymmetric, and therefore we emphasize again that quoting the difference in angle between best fit point and true direction as the angular resolution is somewhat misleading.

Top: simulated best-fit points using the obtained resolution for quality cut (H) of the uncertainties in signal direction, viewing angle and polarization for a point source located at zenith of 50\(^\circ \). For details refer to the text. Bottom: same as the top figure, but for a source located at zenith of 70\(^\circ \). Clearly seen is the zenith dependent shape of the point spread function

Full size image

The orientation of the polarization direction depends on the signal direction, which is different from the neutrino direction, resulting in a ‘bow-tie’-shaped PSF for a neutrino point source, as explained in Sect. 2.5. The geometry of the radio cone, as well as the Vpol-based trigger, limit the allowed orientations of the polarization.

In order to estimate the shape of the PSF, ideally one would want to know the full 3-dimensional probability distribution in signal direction, viewing angle and polarization. Alternatively, one could simulate the distribution of reconstructed events for a given fixed source position. As either option would require a dedicated effort of significant computational cost, we instead proceed as follows. First, we divide the distribution into two subsets based on the quality of the signal direction: ‘well reconstructed’ (\(\varDelta < 2^\circ \)) and ‘poorly reconstructed’ (\(\varDelta > 2^\circ \)). The resolutions for either subset are then assumed to be independent between signal direction, viewing angle and polarization, which we know to be a good approximation for individual events. The PSF is then obtained by folding the resulting resolutions with the distribution of allowed true signal directions for a specific local zenith. We expect the resulting PSF shape to qualitatively agree with the shape that would have been obtained with a dedicated study, and sufficiently accurate to highlight the features of interest.

The PSF for two hypothetical sources, at \(50^{\circ }\) and \(70^{\circ }\) local zenith, are shown in Fig. 11. Only events (hadronic and electromagnetic) passing quality cut (H) are included. For larger incoming zenith angles of the neutrino, a broader part of the cone becomes observable and thereby a larger range of polarization angles. Consequently, the PSF becomes broader and the forbidden region shrinks. Note that these are ‘instantaneous’ PSFs. Except at the poles, a source at fixed declination and right ascension will in general appear at different local zeniths, leading to a slightly larger and more ‘smeared’ time-integrated PSF. However, as long as the detection times are available, using the instantaneous PSF is more accurate and powerful in a source search.

The area of the PSF containing N% of events can be estimated by ‘pixel counting’. As the above procedure results in a space angle distribution that is slightly different from the true distribution, we correct for this by multiplying the area by the ratio \((\varDelta \psi _\textrm{true}/\varDelta \psi _\textrm{PSF})^2\), where \(\varDelta \psi \) denotes the space angle, at each quantile. The resulting areas lie somewhere in between the area expected for a single event and the area of a symmetric PSF, and additionally depend on the local zenith, as demonstrated. For a source at a zenith of 70\(^\circ \) (Fig. 11 bottom, all events (hadronic + electromagnetic) passing quality cut (H)) the median space angle is 2.5\(^\circ \), which would imply a median area of 19\(^\circ \) for a symmetric PSF. The actual PSF area obtained is 14 \(\hbox {deg}^2\) at 70\(^\circ \). At a zenith of 50\(^\circ \), there are two competing effects – the PSF is flatter, but the reconstruction quality is lower (see Fig. 6 and the discussion there). This results, overall, in a larger area than at 70\(^\circ \), albeit with large uncertainties due to the low statistics of events passing quality cut (H) at 50\(^\circ \).

4.4 Systematic uncertainties

To demonstrate the effects of an imperfect knowledge of the detector and the ice target volume for the outlined reconstruction method, we compute an event contour for the same example event shown in Fig. 9, i.e. a typical event remaining after analysis cuts. We take into account three systematic uncertainties, i.e. (1) antenna positions, (2) antenna response and (3) the ice refractive index model. We include an offset in the position of the antenna of 3 cm in any direction, and uncertainty of 10% in the antenna vector effective length and a 10 MHz uncertainty in the frequency behavior of the vector effective length (i.e. a shift in the resonance frequency and the phase behavior). We use Gaussian uncertainties for each individual antenna. For the ice model, we let the parameters \(\varDelta n\) and \(z_{0}\) of the exponential refractive index profile \(n(z) = n - \varDelta n\cdot e^{\frac{z}{z_{0}}}\) vary by \(10\%\), as the deep refractive index can be established very accurately. While further work is needed to correctly treat the systematic uncertainties, we believe that these assumptions provide a reasonable estimate of the current understanding of radio neutrino detectors.

The main impact of the systematic uncertainties for this example event is a worsening of the polarization resolution, from \(\sigma _{68\%}\) = 9.2\(^\circ \) to \(\sigma _{68\%}\) = 11.5\(^\circ \). Additionally, we obtain a bias of 4\(^\circ \) in the reconstructed polarization. This is due to the fact that at low SNR, the arrival time at the Hpol antenna is determined from the arrival times in the nearby Vpol antennas. Uncertainties in the ice model or antenna positions result in a small timing offset between the reconstructed and true pulses, and therefore a systematic underestimation of the Hpol contribution. We emphasize the challenge of the low Hpol signals and state the importance of the development of algorithms to better identify low signal amplitudes in the thermal noise in future work.

The resulting event contour is shown in Fig. 12 resulting in an increase from 13 to 35 \(\hbox {deg}^2\) for the example event. The impact of the mismatch of the model is observed in broader \(\chi ^2_{\text {min}} - \varDelta \chi ^2\) distribution (inset Fig. 12 bottom).

Same event as Fig. 9, including systematic uncertainties

Full size image

4.5 Implications for instrument capabilities

A good pointing resolution for the detected signals is important for neutrino identification, the diffuse flux discovery, cross-section measurements, and source identification. In this section we focus on the latter to illustrate the effect of analysis cuts and angular resolution.

Neutrino sources can be identified in two ways; either through a (temporal and spatial) coincidence with the detection of a photon in the electromagnetic spectrum or a gravitational wave, e.g. [5, 6], or by the detection of an excess of neutrinos that is large enough to identify it above the expected background, e.g. [7]. There are varying analysis strategies for temporal coincidences, where for instance in a very short coincidence window, also a limited angular resolution can lead to a significant detection. However, it is in general understood that a smaller angular uncertainty will always lead to a better chance of finding true coincidences in multi-messenger searches and a reliable arrival direction reconstruction is a prerequisite to a successful triggering of follow-up observations at other wavelengths.

For illustration, we focus here on a potential point source excess, since it is most strongly affected by the angular resolution. A good angular resolution reduces the expected number of background events within the uncertainty region (PSF) of the neutrino (source) and therefore reduces the flux required for source identification.

Identifying point sources by detecting an excess of events above the expected background and implications for the large scale in-ice radio array of IceCube-Gen2 has been extensively studied in [17]. We note that for smaller radio experiments such as RNO-G we rely on very extreme flux models from point sources or nearby sources to be able to identify the source. Nonetheless, we study the effect of the different analysis cuts introduced in Sect. 4, to establish the best search strategy.

We use the software framework toise [58], which allows for the fast calculation of upper limits and discovery potentials of neutrino telescopes based on parameterizations of the detector response in terms of effective volume, angular resolution and energy resolution. The fast performance is due to the usage of an Asimov dataset, meaning that the observed event rates are replaced by the exact mean [59]. Additionally, the test statistic of the null hypothesis (background only) is assumed to be \(\chi ^{2}\)-distributed, i.e. Wilk’s theorem holds [60]. The behavior of the test statistic for the background only scenario can be computed by creating multiple pseudo-experiments using the observed background rates. Background events for UHE neutrinos have not been observed and are as uncertain as the UHE neutrino flux itself, and therefore this approach is as realistic as possible at this stage.

As background for identifying point sources of UHE neutrinos we include three kind of events; (1) the tail of the high-energy neutrino flux, (2) cosmogenic neutrinos, and (3) UHE muons. We include for (1) the diffuse neutrino flux as observed by IceCube, assuming a power-law extension in neutrino energies without a cut-off [61], (2) a cosmogenic neutrino flux resulting from cosmic rays interacting with photon fields assuming a 10\(\%\) proton fraction for the cosmic-ray spectrum observed by the Pierre Auger Observatory [56], and (3) in-ice penetrating UHE muons stemming from cosmic-ray air showers inducing electromagnetic showers while radiating [62]. Each of these contributions are rare and likely on the same level of the expected number of UHE neutrinos, but uncertainties are large. For (1) the flux is measured at low energies (PeV range) with a set of spectral uncertainties, which results in large variations for the EeV regime even if the power-law assumption without cut-off holds. The uncertainties for (2) are dominated by the uncertain UHE cosmic-ray mass composition which influences the neutrino rate prediction by orders of magnitude. Finally for (3) the number of predicted UHE muons in air showers varies due to large uncertainties on the hadronic interaction models, in particular at the highest energies. Therefore, we use the above mentioned background rate as our benchmark background and show results for varying background levels.

Other potential backgrounds for UHE radio detectors are currently under study and are not included in our background estimates. These backgrounds include radio emission from air showers penetrating into the ice and not fully developed in-ice penetrating air showers inducing electromagnetic dense shower cores. Both backgrounds, like the muons, can be reduced if the corresponding air shower can be detected, for example with a shallow component of the hybrid station. We, however, do not include such a cosmic-ray veto in the background assumption.

The impact of the different event groups introduced in Sect. 4 for different background levels is studied in terms of the \(5 \sigma \) discovery potential, i.e. the source flux required for source identification, for a steady point source (10 years) at 20\(^\circ \) declination assuming a neutrino flux of \(E^{-2}\).

Results are shown based on the angular resolution obtained in this work, an energy resolution based on [38], the effective volume based on the specifics of RNO-G in Greenland, using a station effective volume for an RNO-G station and an array of 35 stations [15]. 20\(^\circ \) declination is directly in the field of view of a radio detector based in Greenland.

The point source study modeled with toise is searching 20\(^\circ \) around a source. The framework assumes 1D symmetry for the resolution of a source, i.e. the number of background events is calculated in a circle with radius of the angular resolution. Therefore, the angular resolution required is the PSF, i.e. how a point source is perceived given the imperfect detector. We convert the PSF, obtained as explained in Sect. 4.3, into the area-equivalent symmetric resolution for a source located at a local zenith of 70\(^\circ \). For a station based in Greenland a source at a declination of 20\(^\circ \) has an average local zenith of 70\(^\circ \), which is therefore a good approximation for the size of the PSF.

Fractional improvement of the \(5 \sigma \) discovery potential for different analysis cuts compared to all events with the obtained angular information for a declination of 20\(^\circ \). Note that an event group has a better performance than the benchmark set when the relative discovery potential is lower. We show the analysis cut and two quality cuts (V and H), as discussed in Sect. 4.2. Furthermore, we show the relative discovery potential for the full set using only the shower-maximum resolution. Details are given in the text

Full size image

Our results are visualized in Fig. 13. We show the ratio of the discovery potential for different event groups, improving in angular resolution while reducing the effective volume with their corresponding angular resolution compared to the benchmark event set (full set with obtained resolution). In addition to the event groups introduced in Sect. 4, we also show the results for the full event set with only the shower maximum reconstructed, i.e. the resolution obtained when only exploiting the signal direction. We include this to show the improvement provided through this work.

We conclude that the obtained neutrino direction resolution results in a factor 2–3 improvement in discovery potential compared to reconstructing the shower maximum only. At very low background fluxes, the discovery potential depends only mildly on the number of background events, and thus depends more strongly on the effective volume than the resolution of each event set. As the background flux increases, the reduction in background allowed by restricting to the analysis cut, or eventually one of the quality cuts, becomes favorable, compared to including all events in the analysis. Note, again, that the discovery potential is shown relative to that of all events; for lower background levels or for transient searches the overall discovery potential improves.

In this article, we have described a full reconstruction algorithm for deep in-ice radio neutrino detectors based on the forward-folding principle. We have shown how the signal direction, viewing angle, and polarization combine to uniquely define the neutrino direction. The performance of the algorithm was quantified for each of these parameters, and it was demonstrated that the much larger resolution for polarization leads to strongly asymmetric uncertainty contours for single events. Furthermore, as the alignment of different event contours depends on the interaction vertex of the neutrino rather than its direction, the PSF for a neutrino source in turn exhibits a (zenith-dependent) ’bow-tie’ shape, rather than resembling that of a single event contour. We therefore emphasize that the 1D space angle distribution is not sufficient to adequately describe the reconstruction resolution.

The discovery potential for a point source (for a 20\(^\circ \) search) with this algorithm is shown to improve by a factor \(\approx \) 2 compared to reconstructing the vertex position only. The resolution obtained is shown to depend on the signal strength in the different antennas. For hadronic events, a large subset of events (\(\sim 60\%\)) can be selected in order to obtain a median resolution of \(4.9^\circ \) (space angle) or \(17\ \textrm{deg}^2\) (\(\approx 2.4^\circ \) 1D equivalent). With stricter cuts, a smaller median resolution is obtained at the cost of effective volume. In terms of discovery potential, however, we demonstrate that the loss of statistics generally outweighs the potential improvement in resolution, even for relatively optimistic flux assumptions.

The reconstruction algorithm described in this article is included as part of the open-source NuRadioMC framework, enabling end-to-end simulation and performance studies and optimization for deep in-ice radio neutrino detectors. In future work, improvements are to be expected by improving the identification of pulses with a small signal-to-noise ratio, for which the current strategy of cross-correlation occasionally leads to a failure to correctly identify the pulse position or even to a misreconstruction of the vertex position. In addition, although the performance of the reconstruction for electromagnetic showers is already similar in some aspects (e.g. polarization), further gains are to be expected from a dedicated algorithm for this type of events.