Search for active-sterile antineutrino mixing using neutral-current interactions with the NOvA experiment

This Letter reports results from the first long-baseline search for sterile antineutrinos mixing in an accelerator-based antineutrino-dominated beam. The rate of neutral-current interactions in the two NOvA detectors, at distances of 1 km and 810 km from the beam source, is analyzed using an exposure of $12.51\times10^{20}$ protons-on-target from the NuMI beam at Fermilab running in antineutrino mode. A total of $121$ of neutral-current candidates are observed at the Far Detector, compared to a prediction of $122\pm11$(stat.)$\pm15$(syst.) assuming mixing between three active flavors. No evidence for $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{s}$ oscillation is observed. Interpreting this result within a 3+1 model, constraints are placed on the mixing angles ${\theta}_{24}<25^{\circ}$ and ${\theta}_{34}<32^{\circ}$ at the 90% C.L. for $0.05$eV$^{2} \leq \Delta m^{2}_{41} \leq 0.5$eV$^{2}$, the range of mass splittings that produces no significant oscillations at the Near Detector. These are the first 3+1 confidence limits set using long-baseline accelerator antineutrinos.

This Letter reports results from the first long-baseline search for sterile antineutrinos mixing in an accelerator-based antineutrino-dominated beam. The rate of neutral-current interactions in the two NOvA detectors, at distances of 1 km and 810 km from the beam source, is analyzed using an exposure of 12.51 × 10 20 protons-on-target from the NuMI beam at Fermilab running in antineutrino mode. A total of 121 of neutral-current candidates are observed at the Far Detector, compared to a prediction of 122 ± 11(stat.) ± 15(syst.) assuming mixing only between three active flavors. No evidence forνµ →νs oscillation is observed. Interpreting this result within a 3+1 model, constraints are placed on the mixing angles θ24 < 25 • and θ34 < 32 • at the 90% C.L. for 0.05 eV 2 ≤ ∆m 2 41 ≤ 0.5 eV 2 , the range of mass splittings that produces no significant oscillations at the Near Detector. These are the first 3+1 confidence limits set using long-baseline accelerator antineutrinos.
Studies of neutrinos and antineutrinos from a variety of sources, including accelerators, the atmosphere, the Sun and nuclear reactors [1][2][3][4][5][6][7][8][9][10][11], provide substantial evidence for mixing between the three known active flavors ν τ . However, the observation of antineutrinos in short-baseline oscillations [12,13], reactor [14,15] and gallium-based radiochemical experiments [16,17] found inconsistency with this 3-flavor model. These anomalous results may be explained by extending the mixing model to include more massive neutrino and antineutrino states in addition to (−)  ν 3 with new mass splittings large enough to provide oscillations over much shorter baselines than is possible with the known 3-flavor mass differences. However, observations of the width of the Z 0 -boson resonance at the LEP experiments put the constraint on the number of light antineutrinos participating in the weak interaction to be three [18]. Any additional antineutrinos of mass less than approximately half the Z 0 mass must therefore be sterile [19].
The simplest extension to this 3-flavor model is referred to as the 3+1 model and introduces a single new mass state, (−) ν 4 , with a corresponding sterile flavor state, [20,21]. The matrix can be parametrized as U = R 34 S 24 S 14 R 23 S 13 R 12 [22], where R ij represents a rotation through mixing angle θ ij and S ij represents a complex rotation through angle θ ij and phase δ ij . In addition to the 3-flavor mixing parameters and a new independent mass splitting ∆m 2 41 = m 2 4 − m 2 1 , this model includes three new mixing angles, θ 14 , θ 24 and θ 34 , and two additional phases, δ 14 and δ 24 , which may violate CP symmetry. A number of atmospheric, long-baseline accelerator and reactor neutrino experiments have searched for oscillations outside of the 3-flavor mixing framework and found no evidence of such processes [23][24][25][26][27][28][29], placing constraints on the allowed values of the parameters governing 3+1 oscillations. Equivalent constraints have not yet been published from similar investigations of antineutrinos, the focus of the analysis described in this Letter. This is pertinent as the discrepancies observed in LSND [12] and MiniBooNE [30] involved measurements of oscillations from antineutrino beams, with different behavior reported in subsequent analyses of oscillations of neutrinos [31]. Further investigation is required to understand the oscillations of antineutrinos into potential sterile states.
NOvA can search for evidence of active to sterile oscillations through an analysis of the rate of neutral-current (NC) interactions in its Near Detector (ND) and Far Detector (FD) [29]. When oscillations occur between only the three active flavors, the rate of NC interactions is independent of flavor and is thus unaffected by oscillations. If the antineutrinos transition into an unseen sterile fla-vor state during propagation, a reduction in the number of NC interactions would be observed with a probability given approximately by [32] 1 − P i . To first order, A = sin 2 θ 34 sin 2 2θ 23 and B = 1 2 sin δ 24 sin θ 24 sin 2θ 34 sin 2θ 23 . This approximation assumes ∆ 21 = 0, sin 2 2θ 23 is near-maximal, and considers the fast oscillation regime with limited detector resolution such that sin 2 ∆ 41 = 1 2 ; an exact formalism, including matter effects, is used in the analysis. Studying the disappearance of NC events over a long baseline is therefore sensitive to the 3+1 parameters θ 24 , θ 34 , δ 24 and ∆m 2 41 . Differences between neutrino and antineutrino oscillations due to CP-violation in this approximation are described by the δ 24 phase, and therefore the sign of the B term, and are at most around 10% in the peak antineutrino flux region (around 2 GeV), with a maximum of 20% at higher energies. While an additional asymmetry results from neutrinos gaining an extra effective mass due to forward scattering in matter, this effect is significantly smaller than from the phase.
The ∆m 2 41 parameter determines the frequency of the oscillations. For sufficiently small values, the oscillations develop over longer baselines than the ND and are more rapid than can be resolved at the FD, resulting in an average reduction in the NC rate over the energy spectrum. As the mass splitting increases, oscillations develop over shorter baselines and are evident at the ND. This analysis considers oscillations driven by 0.05 eV 2 ≤ ∆m 2 41 ≤ 0.5 eV 2 , where the NC rate is unaffected at the ND and reduced at the FD.
The NOvA experiment consists of two functionally identical detectors placed in the NuMI beam [33]. The ND is located at Fermilab in Batavia, Illinois, 1 km from the NuMI target and 100 m underground; the FD is located near Ash River, Minnesota, 810 km from the beam source with a 3 m water-equivalent overburden of concrete and barite. The detectors are placed 14.6 mrad off the axis of the beam to produce a narrow neutrino or antineutrino energy distribution peaked at 2.0 GeV at the FD with a width of 0.4 GeV and a high energy tail of around 10% of the peak.
The NuMI beam is produced from 120 GeV protons incident on a mostly carbon target. Two magnetic focusing horns pulsed at 200 kA focus the secondary pions and kaons, which subsequently decay in flight over 705 m, including a 675 m decay pipe, to mainly muons and muon (anti)neutrinos. By selecting the current polarity in the horns, the beam can be run in a mode dominated by either ν µ orν µ . The beam is extracted for 10 µs approximately every 1.3 seconds. A previous anal-ysis searched for evidence of NC-disappearance in a νdominated beam [29]; this study is the first using the NOvAν-dominated dataset. Across the 0-20 GeV energy range, beam interactions at the ND are predicted to consist of 63%ν µ , 35% ν µ and 2% ν e +ν e . In the 0-5 GeV peak energy range, the region which provides greatest sensitivity to θ 34 in this study, beam interactions are 81%ν µ . The analysis signal is NC events, with the sample dominated byν but with some ν component that is treated the same; all charged-current (CC) events are considered background.
The detectors are constructed of highly reflective PVC cells [34] with cross-section 3.9 cm-by-6.6 cm and 15.5 (3.8) m in length for the FD (ND). The cells are formed into 896 (214) planes in the FD (ND), alternating horizontal and vertical relative to the beam axis to enable three-dimensional tracking. Each cell is filled with mineral-oil based liquid scintillator doped with 5% pseudocumene [35] and instrumented with a loop of wavelength-shifting fiber, read out at one end using a 32 pixel avalanche photodiode [36]. The FD and ND have masses of 14 kt and 193 t, respectively. Custom readout electronics shape and digitize the signal. Pulse height and timing of all energy deposits above a given noise threshold are collected in a 550 µs window encompassing the NuMI beam spill at both detectors. Data are additionally collected at the FD using a 10 Hz minimum bias trigger to sample the cosmogenic background.
The analysis uses data recorded between June 2016 and February 2019, corresponding to an exposure of 12.51 × 10 20 protons-on-target (POT) for the FD. The amount of ND data used corresponds to an exposure of 3.10 × 10 20 POT.
A full simulation of the beam, including all materials and interactions in the beamline, is used to predict the flux at the NOvA detectors. Simulation of the antineutrino production from primary proton-induced hadronic showers and subsequent meson transport and decay is provided by a custom made simulation based on GEANT4 [37,38] with a full geometry of the beamline. The flux is then corrected using a version of the Package to Predict the Flux (PPFX) [39] for the NOvA off-axis beam [40], which modifies the particle production using external hadron production data.
Neutrino interactions in the NOvA detectors are simulated by GENIE [41], with modifications performed by NOvA [42]. The particles are propagated through the detector geometry using GEANT4; the charge deposits are converted to scintillation and transported, along with Cherenkov light, to the front-end electronics using a custom made simulation process. The response of the photodetectors is well understood and used to simulate charge collected by the front-end electronics.
NC events are characterized by the lack of a charged lepton exiting the vertex and typically consist of diffuse hadronic activity resulting from energy and momentum transfer from the antineutrino to the nucleus. Interactions are identified by requiring a level of activity consis-tent with these topologies and rejecting CC events with obvious antimuon tracks or positron showers. The NOvA detectors are well-suited to identify NC interactions since the fine granularity results in the hadronic interaction length of 38 cm extending to around 10 cells. This enables a good distinction between π 0 s and positrons, important in selecting NC events containing hadronic activity and rejectingν e CC interactions [43].
The selection was developed on simulated interactions and tested on the ND data, with identical requirements applied to each. The FD has far less overburden than the ND and is exposed to a much higher rate of cosmicinduced interactions; supplementary criteria are therefore applied to FD data to remove the cosmogenic backgrounds. The selection applied in this study, outlined below, is based on the previous neutrino analysis described in [29].
Hits, cells with activity above threshold, are grouped based on their proximity in space and time to form neutrino events candidates. An event vertex is then reconstructed from these hit collections by projecting back to a common start point. Events are required to have a vertex and be fully contained in the detector. This containment requirement stipulates no activity within an outer volume defined for the ND as 25 cm from all sides and for the FD as 100 cm from the top of the detector, 10 cm from the bottom and 50 cm from each of the other sides. Backgrounds arising from neutrino interactions in the material surrounding the ND and resulting in particles, such as neutrons, entering and interacting in the detector are removed by requiring the reconstructed vertex be greater than 100 cm from the side walls and between 150 cm and 1000 cm from the upstream face of the detector.
The primary event classification is provided by a NOvA-custom Convolutional Neural Network (CNN) [44] which identifies interaction topologies from images of the two orthogonal detector views. It is trained separately on simulated FD neutrino-and antineutrino-dominated beam events and cosmic-ray data, and distinguishes the final-state topologies between ν τ CC, NC, and cosmic-induced events. The network learns to extract the characteristic features of the events described above, for example, identifying interactions with charged leptons exiting the vertex as CC, and does not distinguish between interactions from neutrinos and antineutrinos. The separation provided by the CNN for NC events is shown in Fig. 1, and a cut at 0.23 is applied to select events in both detectors. The peak of background events close to 1 is comprised primarily of high energy deep inelastic scattering interactions, or low energy events with minimal detector activity, neither with an obvious lepton to identify. Interactions occurring in the surrounding rock and subsequently entering the ND are removed by requiring the fractional component of the event momentum which is transverse to the beam direction to be less than 80%. At the FD, similar criteria requiring this to be less than 40% for showers with a vertex within 2.4 m from the top of the detector and less than 20% for a vertex inside 2 m are used to remove cosmic-induced interactions. For the FD, candidate antineutrino events must not have significant additional activity within 5 m and 50 ns of the final state particles. Events with activity within 200 cm from the end of the detector must have a sufficient number of downstream planes with hits, compared to the front planes, to ensure showering particles are in the forward direction. A boosted decision tree that considers multiple shower properties is used to further remove cosmic-induced interactions at the FD [45]. Following the application of the full selection, the cosmic background is reduced by over six orders of magnitude. The selected sample in the FD simulation is 78% pure. Two-thirds of the remaining backgrounds result from CC interactions where most of the energy is transferred to the hadronic system, and the rest originate from cosmicinduced interactions. The NC-candidate events selected from beam neutrino interactions are 79%ν and 21% ν in the 0-5 GeV peak flux region. The antineutrino-candidate energy is estimated for NC events by summing the calibrated energy deposited in the interaction and applying a linear scaling from simulation to correct for the unseen outgoing energy. The final reconstructed energy has a resolution of around 30%.
The spectrum of selected NC-candidate events at the ND is shown in Fig. 2. From this, the predicted spectrum at the FD is produced following a data-driven approach in which the ND data are used to constrain the selected simulated events before propagating the resulting distribution to the FD [46]. This procedure uses the simulated ratio between the spectra at the two detectors and takes into account geometric differences and effects of the beam dispersion. The flavor composition and energy of the corrected ND NC-candidate spectrum is inferred from simulation and used to apply oscillations to the events. This approach is effective in partially cancelling any correlated uncertainties between the two detectors but results in sensitivity only to active-to-sterile oscillations which develop over larger distances than the ND baseline. In a 3+1 model these slower oscillations are driven by mass splittings in the range 0.05 eV 2 ≤ ∆m 2 41 ≤ 0.5 eV 2 . The dominant uncertainties are summarized in Table I. The largest are attributed to the calibration process which determines the energy scale of the charge collected in each detector cell. An absolute calibration uncertainty of 5% is used, consistent with the agreement with simulation of proton response in quasielastic-like ND CC interactions. A calibration shape uncertainty to cover charge deposited near the poorly-modeled detector edges is evaluated by considering the differences between calibrated and simulated deposits. These uncertainties are applied both as correlated absolute offsets at both detectors and separately as relative differences between them, since they are determined using just the ND. Reducing this large calibration uncertainty is a priority for subsequent analyses and there are various ongoing efforts, such as improving the analysis procedure and utilizing understanding gained from a test beam program [47].
The simulation contains a modeling of the detector response, including the creation of scintillation and Cherenkov light, its transport to the front-end electronics and subsequent processing. An uncertainty of 10% in the normalization of the scintillation light is determined from differences between data and simulation in the number of photoelectrons collected by through-going minimally-ionizing particles. The residual uncertainties are quantified using alternate predictions where the scintillation and Cherenkov photon yields in the model are altered within the tolerance of agreement with the ND data while holding the muon response fixed, since that is set by the calibration procedure.
Neutrino interaction model uncertainties are estimated using the GENIE reweighting framework, with additional uncertainties on the modifications of the various interaction models performed by NOvA [42]. An additional uncertainty is considered to account for the differences between the selected spectrum with and without applying these NOvA reweighting modifications. A further 60% uncertainty on the ν τ CC cross-section is allowed, consistent with the most recent measurements [48][49][50][51].
Uncertainties in particle propagation when simulating the beam are provided by the PPFX reweighting framework [39]. Additional transport effects are included by considering the uncertainties in the positioning and modelling of all beamline components including target, focusing horns and decay volume. The high-energy antineutrinos in the spectrum are produced from kaons, rather than pions, and these hadrons have a 30% uncertainty on their simulated production. The magnitude of this is motivated by the observed discrepancy with the PPFX prediction in previous MINOS+ studies [23].
The fraction of energy carried by neutrons that results in visible signals is uncertain. This is quantified by considering the difference in reconstructed energy between data and simulation for quasielastic-like CC events with one µ + and a single additional reconstructed object in the final state. This second particle is selected using criteria consistent with having a neutron parent based on energy, displacement from the vertex and angle. The observed deficit in reconstructed energy for these neutrons is used to estimate the missing energy scale. An uncertainty on the NC-candidate event spectrum is then evaluated using this information to determine the underestimation of the reconstructed antineutrino energy from neutrons.
Each systematic uncertainty is applied to the simulation and used to produce oscillated FD spectra using the data-driven approach previously described. The resulting changes in the event spectrum are taken as the residual energy shifts and included for both signal and background events as additional parameters in the fit. The total uncertainty on simulated FD events following this procedure is shown in Fig. 3. The event selection and treatment of systematic uncertainties were finalized based on simulation before comparing with the FD data. Upon analyzing the data, 121 NC-candidate events were observed, consistent with ex-   Table II. World-average 3-flavor oscillation parameters were used [18], assuming the most conservative case of normal mass ordering and upper θ 23 octant. The observed spectrum in Fig. 3 shows good agreement between data and prediction across the range of energies.
To provide a model-independent metric of agreement with an expectation assuming oscillations between only the three active antineutrino species, the ratio R NC ≡ N − P bkg P NC [52] is used to compare the number of data events, N , with the predicted number of signal and background events, P NC and P bkg respectively. This is equal to unity in the case of no NC-disappearance. The obtained value of 0.99 ± 0.12(stat.) +0.14 −0.16 (syst.) is consistent with no oscillations involving a sterile antineutrino.
A fit considering rate and shape is performed on the spectrum in Fig. 3 to limit the allowed values of the parameters governing sterile oscillations in the 3+1 framework. The metric is minimized, varying the number of events, P i , predicted  [25] and IceCube DeepCore [27]. The limits from Super-Kamiokande and IceCube are set using atmospheric neutrinos, a sample which is majority neutrinos but with a non-negligible antineutrino component. The NOvA limits are the first set from studying an antineutrino-dominated sample.
in each FD bin, i, and constraining the agreement with data, N i , using Gaussian-distributed penalty terms, U j , for each systematic uncertainty j, each with width σ j . The mixing angle θ 14 and the corresponding phase δ 14 are set to zero, consistent with results from solar and reactor experiments [24] and unitarity considerations [53], and the δ 24 CP-violating phase was included as a nuisance parameter in the fit, allowed to float freely without penalty. The δ 24 parameter weakens the limit for large values of both θ 24 and θ 34 but does not significantly impact the individual limits. The fit is performed at ∆m 2 41 = 0.5 eV 2 but the result is valid for the full range 0.05 eV 2 ≤ ∆m 2 41 ≤ 0.5 eV 2 . Confidence limits are computed using the unified approach of Feldman-Cousins [54] and are shown in Fig. 4.
Profiling over each sterile mixing angle in turn provides limits of θ 24 < 25 • and θ 34 < 32 • at 90% C.L. Table III shows these constraints, together with the inferred ma-trix element values, alongside results from other experiments. This result is compatible with world limits from analyses reporting no evidence of sterile neutrino oscillations, shown in Fig. 4, and is the first from an analysis of antineutrinos over a long-baseline.
To conclude, 121 neutral-current candidates are observed at the Far Detector compared with 122 ± 11(stat.) ± 15(syst.) events predicted assuming mixing only occurs between active species. This is the first accelerator-based search for oscillations into sterile neutrinos over long-baselines using an antineutrinodominated beam, and no evidence of a depletion of the active flavors is observed. The measurement reported in this Letter enhances the understanding of antineutrino oscillations in the context of sterile neutrino anomalies. Looking forward, NOvA plans to increase its antineutrino dataset by a factor of around 2.5, complementing a comparable neutrino sample, which will facilitate improved sterile neutrino searches and enable more precise measurements of the governing parameters.  [23], T2K [28], Super-Kamiokande [25], IceCube [26] and IceCube-DeepCore [27]. The limits are shown for ∆m 2  grateful for the contributions of the staffs of the University of Minnesota at the Ash River Laboratory, and of Fermilab.