Ultrahighenergy neutrinos and cosmic rays are excellent probes of astroparticle physics phenomena. For astroparticle physics analyses robust and accurate reconstruction of signal parameters like arrival direction and energy is essential. Current reconstruction methods ignore bintobin noise correlations which limits reconstruction resolution and so far has prevented calculations of eventbyevent uncertainties. In this work we present a likelihood description of neutrino or cosmicray signals in a radio detector with correlated noise as present in all neutrino and cosmicray radio detectors. We demonstrate with a toymodel reconstruction that signal parameters such as energy and direction including eventbyevent uncertainties with correct coverage can be obtained. Additionally by correctly accounting for correlations the likelihood description constrains the bestfit parameters better than alternative methods and thus improves experimental reconstruction capabilities.
['Astrophysics - Instrumentation and Methods for Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Likelihood reconstruction of radio signals of neutrinos and cosmic rays
{'No Header': '<!-- image -->  \n<!-- image -->', 'Martin Langgård Ravn, 𝑎, ∗ Chistrian Glaser, 𝑎 Thorsten Glüsenkamp 𝑎 and Alan Coleman 𝑎': '𝑎 Dept. of Physics and Astronomy, Uppsala University, Box 516, S-75120 Uppsala, Sweden  \nE-mail:  \nmartin.ravn@physics.uu.se  \nUltra-high-energy neutrinos and cosmic rays are excellent probes of astroparticle physics phenomena. For astroparticle physics analyses, robust and accurate reconstruction of signal parameters like arrival direction and energy is essential. Current reconstruction methods ignore bin-to-bin noise correlations, which limits reconstruction resolution and so far has prevented calculations of event-by-event uncertainties. In this work, we present a likelihood description of neutrino or cosmic-ray signals in a radio detector with correlated noise, as present in all neutrino and cosmic-ray radio detectors. We demonstrate with a toy-model reconstruction that signal parameters such as energy and direction, including event-by-event uncertainties with correct coverage, can be obtained. Additionally, by correctly accounting for correlations, the likelihood description constrains the best-fit parameters better than alternative methods and thus improves experimental reconstruction capabilities.  \n10th International Workshop on Acoustic and Radio EeV Neutrino Detection Activities  \n11 - 14 June, 2024  \nThe Kavli Institute for Cosmological Physics, Chicago, IL, USA', '1. Introduction': 'In this contribution, we present a probabilistic description of voltage time-traces recorded in radio detectors with band-limited noise and a deterministic neutrino or cosmic-ray signal. The probabilistic description can be viewed as the likelihood for a signal given the measured traces and is thus the optimal objective function for reconstruction. We demonstrate with a toy-model that, when using the likelihood description, the signal parameters can be reconstructed with uncertainties, the uncertainties give correct coverage, and the reconstruction resolution improves compared to previous methods.  \nIn recent years, radio detectors have proven to be powerful tools for probing ultra-high-energy astroparticle phenomena. The radio detection technique has been successfully employed in cosmicray air shower observatories, such as the Pierre Auger Observatory and LOFAR, and has yielded accurate estimation of arrival direction, shower maximum, and energy with precision competitive to other established methods [1]. In the field of neutrino astronomy, radio detectors have emerged as the most promising approach to detect the yet undiscovered ultra-high-energy neutrino flux at PeV energies and above [2]. Several prototype in-ice radio arrays, such as ARIANNA and ARA, along with balloon experiments, like ANITA, have demonstrated the viability of this technique. Furthermore, the construction of larger in-ice arrays, such as RNO-G, is ongoing.  \nIn any astroparticle physics analysis involving radio-detected signals, robust and accurate estimation of the signal parameters is crucial. For cosmic-ray and neutrino signals, the parameters of interest are the energy of the primary particle and its arrival direction. Additionally, for cosmicray studies, the shower maximum is an important proxy for the composition.  \nSeveral successful methods for reconstructing neutrino or cosmic-ray signals in radio detectors have been presented in recent years. These are either based on summary statistics calculated for each antenna observing the signal (e.g. signal arrival time, maximum amplitude, and fluence), or they use the full time-trace measured by each antenna. In general, the time-traces hold all available information and should result in optimal reconstruction, however, since antenna responses and other detector effects can be challenging to model and measure, methods using summary statistics sometimes result in a more stable reconstruction.  \nThe current method for in-ice radio detectors is the forward folding technique [3], in which the radio emission of the neutrino or cosmic-ray interaction is calculated, propagated through the given medium, and folded through the antenna and hardware response to obtain the signal as it would appear in the recorded time-trace. The predicted signal is then compared to a measured time-trace, where bin-by-bin the chi-square between the predicted and measured voltage is calculated and summed over all bins and antennas:  \n𝜒 2 = 𝑛 ant ∑︁ i = 1 𝑛 𝑡 -1 ∑︁ k=0 [ 𝑉 𝑖,𝑘 -𝜇 𝑖,𝑘 ( 𝜃 )] 2 𝜎 2 𝑘 , (1)  \nwhere 𝑉 𝑖,𝑘 is the measured voltage in the 𝑖 th antenna and 𝑘 th time bin, 𝜇 𝑖,𝑘 ( 𝜃 ) is the predicted signal parameterized by 𝜃 , 𝜎 𝑘 is the root mean square (RMS) of the noise in the 𝑘 th antenna, 𝑛 ant is the number of antennas, and 𝑛 𝑡 is the number of time bins. This 𝜒 2 is then used as an objective function to minimize in a reconstruction of signal parameters, which has been shown to work well for reconstruction of electrical fields from cosmic rays and neutrino signals [3-5]. While references  \n[3, 4] are simulation studies, reference [5] used the forward folding technique to reconstruct cosmic rays measured by the ARIANNA detector and demonstrates the applicability of the method under experimental conditions.  \nHowever, the 𝜒 2 expressed in Equation (1) can not be used to estimate uncertainties on the reconstructed parameters on an event-by-event basis. The reason for this is that the noise present in radio detectors is not white noise, but can, in general, be characterized as band-limited, which results in correlations between neighboring bins in the time domain. Since the correlations are not taken into account in the 𝜒 2 , it is not 𝜒 2 -distributed and can thus not be used to estimate uncertainties. Instead, the correct probabilistic description of band-limited noise plus a deterministic signal, which takes the correlations in the noise into account, is needed to calculate the likelihood for a signal. We present such a model in this work.', '2. Probabilistic noise model': "We model the band-limited noise plus a deterministic signal in a time-trace as a multivariate normal distribution with dimensionality equal to the length of the trace, 𝑛 𝑡 . The probability density of measuring any trace 𝒙 = ( 𝑉 0 , ..., 𝑉 𝑛 𝑡 -1 ) is then:  \n𝑝 GLYPH<0> 𝒙 ; 𝝁 ( 𝜃 ) , 𝚺 GLYPH<1> = 1 √︁ ( 2 𝜋 ) 𝑛 𝑡 | 𝚺 | exp GLYPH<18> -1 2 GLYPH<0> 𝒙 -𝝁 ( 𝜃 ) GLYPH<1> T 𝚺 -1 GLYPH<0> 𝒙 -𝝁 ( 𝜃 ) GLYPH<1> GLYPH<19> , (2)  \nwhere 𝝁 ( 𝜃 ) = ( 𝜇 0 ( 𝜃 ) , ..., 𝜇 𝑛 𝑡 -1 ( 𝜃 )) is the signal (zeros if no signal is present) parameterized by 𝜃 , 𝚺 is the covariance matrix of the noise, and | 𝚺 | is its determinant. The covariance matrix can be estimated from many traces consisting purely of noise via the estimate of the autocovariance function:  \nCov ( 𝑡 𝑘 , 𝑡 𝑘 ' ) = 1 𝑁 𝑁 ∑︁ 𝑛 = 1 ( 𝑥 𝑛,𝑘 -¯ 𝑥 𝑘 )( 𝑥 𝑛,𝑘 ' -¯ 𝑥 𝑘 ' ) , (3)  \nwhere 𝑁 is the number of traces and ¯ 𝑥 𝑘 is the mean of the 𝑘 th time bin, which is zero for traces of pure noise. By assuming the covariance matrix is a symmetric circulant matrix, the elements only depend on Δ 𝑡 𝑚 = 𝑡 𝑘 ' -𝑡 𝑘 , where 𝑚 = 𝑘 ' -𝑘 , and we can average over elements with the same Δ 𝑡 𝑚 :  \n⟨ Cov ( Δ 𝑡 𝑚 )⟩ = 1 𝑛 𝑡 𝑛 𝑡 -1 ∑︁ 𝑘 = 0 𝑛 𝑡 -1 ∑︁ 𝑘 ' = 0 𝛿 𝑘, ( 𝑘 ' -𝑚 ) Cov ( 𝑡 𝑘 , 𝑡 𝑘 ' ) . (4)  \nwhere 𝛿 is the Kronecker delta. The elements of the covariance matrix are then:  \n𝚺 𝑘𝑘 ' = ⟨ Cov ( Δ 𝑡 𝑘 ' -𝑘 )⟩ . (5)  \nThecirculant structure of the covariance matrix reflects that the noise correlations are invariant under time translation and are periodic. Although the periodic assumption is not physically accurate, since it implies correlations between the first and last bins of a trace, it has minimal impact on this work as the signals considered can be assumed not to be close to the edges of the traces. Alternatively, a Toeplitz structure can be employed, which does not impose the periodic constraint. However, the circulant structure is more numerically stable and simplifies the implementation, and is thus  \n<!-- image -->  \nFigure 1: Left: Band-limited frequency spectrum, which is typical for radio detectors. Frequency amplitudes below 1 h of the maximum height of the spectrum are set to 0 to avoid numerical instabilities. Right: Covariance matrix corresponding to the spectrum represented as an image (zoom in on the first 50 ns).  \n<!-- image -->  \nused in the following. An example of a typical frequency spectrum for radio detectors and the corresponding covariance matrix estimated from 10000 realizations of noise are shown in Figure 1.  \nWe note that the multivariate normal distribution (Equation (2)) depends on the inverse of the covariance matrix, which can only be calculated if it is full rank. The covariance matrix is, however, not full rank when any frequency amplitudes are exactly zero. In this case, the multivariate normal distribution is calculated by replacing the inverse with the Moore-Penrose pseudo inverse 𝚺 -1 → 𝚺 + , the determinant with the pseudo-determinant | 𝚺 | → | 𝚺 | + , and 𝑛 𝑡 with two times the number of non-zero frequency amplitudes 𝑛 𝑡 → 𝑛 dof , which is the number of degrees of freedom.", '3. Likelihood reconstruction': "The multivariate normal probability of Equation (2) can be viewed as the likelihood for a signal given a trace. When the signal is present in many antennas, the full likelihood for the signal is the product of the likelihood for each antenna in the detector. The minus two log-likelihood is then:  \n-2 ln L( 𝝁 ( 𝜃 ) ; 𝒙 , 𝚺 ) = 𝑛 ant ∑︁ 𝑖 = 1 GLYPH<0> 𝒙 𝑖 -𝝁 𝑖 ( 𝜃 ) GLYPH<1> T 𝚺 -1 𝑖 GLYPH<0> 𝒙 𝑖 -𝝁 𝑖 ( 𝜃 ) GLYPH<1> + 𝑐𝑜𝑛𝑠𝑡. (6)  \nwhere 𝒙 , 𝝁 𝑖 ( 𝜃 ) , and 𝚺 𝑖 are the trace, signal, and covariance matrix of the 𝑖 th antenna, and the constants can be ignored since, in general, only delta log-likelihoods are of interest, i.e., likelihood ratios. To demonstrate that the likelihood description can be used to reconstruct signals and estimate uncertainties on the reconstructed parameters, we perform a toy-model reconstruction of a simplified neutrino signal in a toy-model in-ice radio detector. The neutrino signal calculation and detector simulation are performed using NuRadioMC [6] and NuRadioReco [3].  \nThe toy-model detector layout is a typical 'deep' station layout similar to the design of the ARA detector at the South Pole [7] and RNO-G in Greenland [8], and similar to the deep detector component foreseen for IceCube-Gen2 [9]. It consists of 12 antennas in total, located in three holes in the ice. Each string has both a vertically polarized (VPol) and a horizontally polarized (HPol)  \nantenna at 80 m depth and the same at 100 m depth. The coordinates of the strings in the horizontal plane ( 𝑥 , 𝑦 ) are (0 m, 0 m), (20 m, 0 m), and (0 m, 20 m) which form a right isosceles triangle.  \nThe Askaryan radio emission from the neutrino interaction is calculated using the Alvarez2009 [10] parameterization as implemented in NuRadioMC, assuming a hadronic shower. The most computationally intensive part of the signal calculation is the propagation of the radio signal through the ice. For ease of computation, we thus neglect the refraction of radio waves in non-uniform media and assume the signal travels in a straight line from the interaction vertex to the detector. This simplification does not affect the conclusions of this work. After the signal is propagated through the ice, it is folded through the antenna response, and noise is added to the resulting trace. The noise is generated with the frequency spectrum shown in Figure 1, which is typical for radio detectors. We set the noise temperature to 300 K. The traces are 𝑛 𝑡 = 1024 samples long and have a sampling rate of 1 . 6 GHz.  \nThe resulting neutrino signal observed in the detector depends on seven parameters: the shower energy, 𝐸 shower , the neutrino arrival direction zenith and azimuth angles, 𝜃 𝜈 and 𝜙 𝜈 , the interaction vertex position in spherical coordinates, 𝑟 vertex , 𝜃 vertex and 𝜙 vertex , and the signal arrival time, 𝑡 0 . These are the parameters we aim to estimate in a reconstruction. In this work, the results for three representative example events are shown, labeled event 1, 2, and 3. The Monte Carlo truth values of the parameters for the three events are shown in Table 1. The traces in all 12 antennas for Event 3 with one realization of noise is shown in Figure 2.  \nTable 1: True parameter values for the three example events shown in this work.  \nFigure 2: Traces for the 12 antennas of the toy-model in-ice radio detector with a neutrino signal (Event 3 in Table 1) and band-limited noise (left panel in Figure 1).  \n<!-- image -->  \nThe reconstruction is performed by minimizing the -2 ln L( 𝝁 ( 𝜃 ) ; 𝒙 , 𝚺 ) of Equation (6) with respect to the neutrino signal parameters 𝜃 = ( 𝐸 shower , 𝜃 𝜈 , 𝜙 𝜈 , 𝑟 vertex , 𝜃 vertex , 𝜙 vertex , 𝑡 0 ) using the iminuit [11] MIGRAD algorithm. To investigate the ideal case where the global minimum of the likelihood can be found, the parameters are initialized at the true values in the minimization process.  \nTo estimate uncertainties on a subset of the reconstructed parameters, we perform a profile likelihood scan. In many physics analyses, the primary parameters of interest are the neutrino energy or the arrival direction, while the remaining parameters can be treated as nuisance parameters and profiled over. In this study, we focus specifically on the neutrino arrival direction. The profile likelihood scan procedure is as follows: First, the reconstruction is performed with all parameters free to find the global minimum of the -2 ln L , denoted as -2 ln L( ˆ 𝜃 ) . Then, a grid is defined in 𝜃 𝜈 and 𝜙 𝜈 , and for each point in the grid, a minimization is performed with 𝜃 𝜈 and 𝜙 𝜈 fixed at the corresponding values. This yields a -2 ln L value for each alternative hypothesis, which we denote -2 ln L( 𝜃 ' ) . According to Wilks' theorem [12] the -2 Δ ln L = -2 [ ln L( 𝜃 true ) -ln L( ˆ 𝜃 )] between the best fit and true parameter values follows a 𝜒 2 -distribution with 2 degrees of freedom when profiling over the nuisance parameters. The -2 Δ ln L = -2 [ ln L( 𝜃 ' ) -ln L( ˆ 𝜃 )] between the best fit and the alternative hypothesis can then be used to determine a confidence level for each point in the grid and draw uncertainty contours for the neutrino arrival direction. The reconstructed neutrino arrival directions, along with uncertainty contours derived from the profile likelihood scans for the three example events with one realization of noise each, are displayed in the top row of Figure 3.  \nFigure 3: Top: Profile likelihood scans for the three example events listed in Table 1 with one realization of noise each along with the reconstructed arrival direction, 1 𝜎 uncertainty contours, and the true arrival direction. Bottom: Coverage for 100 trials of the same neutrino signals with different realizations of noise for the three example events using the -2 ln L expressed in Equation (6) presented in this work along with the coverage using the 𝜒 2 with no correlation taken into account shown in Equation (1).  \n<!-- image -->  \nThe validity of the estimated uncertainties can be verified by repeating the reconstruction many times using the same true parameter values and different realizations of noise. For any confidence level, we can then count how many times the true value would have been within the corresponding uncertainty contours, which gives the coverage. The confidence versus coverage is shown in the bottom row of Figure 3 for 100 trials of each of the three example events. The figure shows that the coverage is correct for any confidence level for the true parameter values tested here. Using the likelihood description presented in this work, we can thus reconstruct neutrino signals in radio detectors with band-limited noise and correctly estimate uncertainties on the reconstructed parameter values.  \nIn order to investigate the impact of using the likelihood description, we repeat the reconstruction using the 𝜒 2 as expressed in Equation (1) as the objective function which does not take the correlations between bins into account. In this case, the reconstruction does successfully return best-fit parameters, however, two issues are present. First, if the 𝜒 2 is assumed to be 𝜒 2 -distributed and used to estimate uncertainties on the reconstructed parameters, the coverage is completely wrong, which is shown in the dotted lines in the bottom row of Figure 3. Hence, uncertainties can not be estimated for measured neutrinos, as resimulation with different noise realizations is only possible in simulation studies. Second, the reconstructed parameters for 100 trials have a significantly larger spread around the true values. The spread of the seven reconstructed parameters for 100 trials of the three example events using the two objective functions is shown in Table 2. This result emphasizes that to achieve the best possible reconstruction resolution of signals in radio detectors, the correct likelihood description is needed, which employs the information about the correlations in the noise.  \nTable 2: Spread (standard deviation) of the reconstructed parameter values for 100 trials of the three example events listed in Table 1 minimizing the -2 ln L presented in Equation (2) in this work compared to minimizing the 𝜒 2 expressed in Equation (1), which does not take correlation into account.", '4. Summary and outlook': 'In this work, we have presented a likelihood description of a deterministic signal in a radio detector with band-limited noise. The likelihood allows us to obtain uncertainties on the reconstructed parameters, which was not possible with previous approaches. We have demonstrated with a toy-model study using this likelihood description that the uncertainties have correct coverage. We also find that the Likelihood yields smaller uncertainties for in-ice neutrino reconstruction compared to the previously employed 𝜒 2 function that does not take into account correlations.  \nThe likelihood description is generally a strong tool with several additional applications currently being investigated. Since the likelihood ratio is the best test statistic to compare two hypotheses, it can be used to separate signal from background. This fact can potentially be used to improve the identification of cosmic-ray signals in in-ice radio detectors, like in the cosmic-ray search recently performed with the ARIANNA detector [5]. Second, the likelihood description provides a way to calculate a goodness-of-fit for a reconstructed signal. For a reconstructed signal that describes the measured trace well, the -2 ln L in Equation (2) (ignoring the constants) should be 𝜒 2 -distributed with degrees of freedom equal to the number of samples in the trace from which a p-value can be calculated. Finally, the covariance matrix discussed in Section 2 can be used to calculate the Fisher information matrix for a detector and set of signal parameters. The inverse Fisher information matrix gives an estimate of the uncertainties and correlations of the signal parameters without the need to perform any reconstruction. The likelihood description presented in this work is thus potentially also useful in the prospect of end-to-end detector optimization.', 'References': '- [1] F. G. Schröder, Radio detection of Cosmic-Ray Air Showers and High-Energy Neutrinos , Prog. Part. Nucl. Phys. 93 (2017) 1-68, [ 16/zero.alt37./zero.alt38781 ].\n- [2] S. Barwick and C. Glaser, Radio Detection of High Energy Neutrinos in Ice , The Encyclopedia of Cosmology, chapter 6, pp. 237-302 (2023) , [ 22/zero.alt38./zero.alt34971 ].\n- [3] C. Glaser, A. Nelles, I. Plaisier, C. Welling, S. W. Barwick, D. García-Fernández et al., NuRadioReco: A reconstruction framework for radio neutrino detectors , Eur. Phys. J. C 79 (2019) 464, [ 19/zero.alt33./zero.alt37/zero.alt323 ].\n- [4] I. Plaisier, S. Bouma and A. Nelles, Reconstructing the arrival direction of neutrinos in deep in-ice radio detectors , Eur. Phys. J. C 83 (2023) 443, [ 23/zero.alt32./zero.alt3/zero.alt3/zero.alt354 ].\n- [5] A/r.pc/i.pc/a.pc/n.pc/n.pc/a.pc collaboration, A. Anker et al., Measuring the polarization reconstruction resolution of the ARIANNA neutrino detector with cosmic rays , JCAP 04 (2022) 022, [ 2112./zero.alt315/zero.alt31 ].\n- [6] C. Glaser et al., NuRadioMC: Simulating the radio emission of neutrinos from interaction to detector , Eur. Phys. J. C 80 (2020) 77, [ 19/zero.alt36./zero.alt3167/zero.alt3 ].\n- [7] RNO-G collaboration, P. Allison et al., Design and Initial Performance of the Askaryan Radio Array Prototype EeV Neutrino Detector at the South Pole , Astropart. Phys. 35 (2012) 457-477, [ 11/zero.alt35.2854 ].\n- [8] RNO-G collaboration, J. A. Aguilar et al., Design and Sensitivity of the Radio Neutrino Observatory in Greenland (RNO-G) , JINST 16 (2021) P03025, [ 2/zero.alt31/zero.alt3.12279 ].\n- [9] I/c.pc/e.pcC/u.pc/b.pc/e.pc-G/e.pc/n.pc2 C/o.pc/l.pc/l.pc/a.pc/b.pc/o.pc/r.pc/a.pc/t.pc/i.pc/o.pc/n.pc collaboration, I.-G. Collaboration, Icecube-gen2 technical design report , 2023.\n- [10] J. Alvarez-Muniz, W. R. Carvalho, Jr., M. Tueros and E. Zas, Coherent Cherenkov radio pulses from hadronic showers up to EeV energies , Astropart. Phys. 35 (2012) 287-299, [ 1/zero.alt3/zero.alt35./zero.alt3552 ].\n- [11] H. Dembinski and P. O. et al., scikit-hep/iminuit , .\n- [12] S. S. Wilks, The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses , Annals Math. Statist. 9 (1938) 60-62.'}