Constraining the NuMI neutrino flux using inverse muon decay reactions in MINERvA

Neutrino oscillation experiments [1, 2, 3, 4] depend on measurements of neutrino interactions at a near detector as a companion measurement that probes the flux and neutrino interaction cross sections that affect the experiment. However, there are significant uncertainties both in the cross sections for neutrino interactions and in the reconstruction of neutrino energies of most reactions observed at near detectors. These uncertainties make it difficult to use only measurements at a near detector to measure the neutrino flux and separate it from the effects of neutrino interactions.

One partial solution to this problem is to measure scattering of neutrinos from atomic electrons. Such scattering is accurately predicted in the Standard Model, with uncertainties of a per cent or less primarily due to hadronic effects in radiative corrections [5]. These reactions then provide a measurement of the flux which is independent of interaction uncertainties and can help to break degeneracies between those interaction uncertainties and uncertainties in predictions of the flux. This technique has been demonstrated by the MINERvA experiment in $\nu e^{-}\to\nu e^{-}$ scattering [6, 7] and has been studied for application in the future DUNE [3] experiment [8].

Another neutrino-electron scattering reaction is inverse muon decay (IMD), $\nu_{\mu}e^{-}\to\nu_{e}\mu^{-}$. The IMD process has a threshold energy of $E_{\mathrm{min}}=\frac{m_{\mu}^{2}-m_{e}^{2}}{2m_{e}}\approx 11$ GeV, and a total cross section given at tree level by [9]

For the purposes of constraining neutrino flux, IMD is only sensitive to a single neutrino type in the beam, muon neutrinos, and is only initiated by neutrinos above the threshold. From just the spectrum of muons alone, there is only a weak correlation between muon energy and neutrino energy. Therefore, the number of IMD events measures some weighted integral of $\nu_{\mu}$ over the reaction threshold. For the NuMI neutrino beam, whose neutrino-dominated (“forward horn current” or FHC) and anti-neutrino-dominated “reverse horn current” or RHC) fluxes are shown in Fig. 1, the focusing peak is below the threshold so IMD is sensitive only to the energies greater than the focusing peak, the “high-energy tail”, of the beam. This tail has a large contribution from neutrinos which are unfocused or under-focused by the beam optics [10, 11].

Backgrounds to the measurement come almost entirely from high-energy neutrino $\nu_{\mu}$ quasielastic scattering on bound neutrons in nuclei, with small contributions from multi-nucleon and inelastic processes. Background models described below will be improved with constraints from “sideband” samples at lower $E_{\mu}$ and higher $\theta_{\mu}$ than the IMD signal.

The MINERvA experiment employs a fine-grained tracking detector for recording neutrino interactions produced by the NuMI beamline at Fermilab [11, 12]. Neutrinos are created by directing 120 GeV protons from the Main Injector onto a graphite target. The resulting charged pions and kaons are focused by two magnetic horns. Choice of the polarity of the current in the magnetic horns gives either the FHC or RHC beams, as defined above, and this analysis uses data from both beams. Approximately 97$\%$ of the muon neutrinos that reach the MINERvA detector are produced by pion decay; the remainder are the result of kaon decay [11, 12]. At the largest neutrino energies, the fraction of neutrinos from kaons increases. For neutrinos produced from the highest energy pions and kaons, the focusing from the horns is generally ineffective, and so the numbers of high-energy $\nu_{\mu}$ in the FHC and RHC beams, particularly those originating from $\pi^{+}$ decays, are similar.

The MINERvA detector [13] consists of 120 hexagonal modules that create an active tracking volume preceded by a set of passive nuclear targets. This result includes only those interactions in the active tracking volume with a fiducial mass of 5.48 tons. The active target volume is surrounded by electromagnetic and hadronic calorimeters. Each tracking module has two planes composed of triangular polystyrene scintillator strips with a 1.7 cm strip-to-strip pitch. For three-dimensional reconstruction, planes are oriented in three different directions, 0^∘ and $\pm$ 60^∘ relative to the vertical axis of the detector. The downstream and side electromagnetic calorimeters consist of alternating layers of scintillator and 2 mm-thick lead planes. The downstream hadronic calorimeter consists of alternating scintillator and 2.54 cm-thick steel planes. Multi-anode photomultiplier tubes read out the scintillator strips via wavelength-shifting fibers. The timing resolution of the readout electronics is 3.0 ns and sufficient to separate multiple interactions within a single NuMI beam spill.

Muons that originate in MINERvA from IMD travel entirely through MINERvA into the MINOS near detector [14] located 2 m downstream of the MINERvA detector. In MINOS, their momentum and electric charge are measured by a magnetized spectrometer composed of scintillator and iron.

This analysis uses data that correspond to $10.61\times 10^{20}$ protons on target (POT) in the FHC configuration and $11.24\times 10^{20}$ POT in the RHC configuration taken between September 2013 and February 2019. The beam focusing configuration and target are that of the “medium energy” beam provided for the NOvA experiment.

A GEANT4-based simulation of the NuMI beamline is used to predict the neutrino flux. To improve the prediction, the simulation is reweighted as a function of pion kinematics to correct for differences between the GEANT4 [15] prediction and hadron production measurements of 158 GeV protons on carbon from the NA49 experiment [16] and other relevant hadron production measurements. A description of this procedure is found in Ref. [12]. The in situ measurement of neutrino scattering off atomic electrons described in Ref. [7], is not used in this analysis to constrain the flux prediction. This measurement and the neutrino-electron elastic scattering measurements give independent constraints which may be combined.

The kinematics of the parent mesons for IMD neutrinos in both FHC and RHC, as predicted by the simulation, are shown in Fig. 2. This simulation predicts two dominant populations of mesons that produce neutrinos with sufficient energy to contribute to the inverse muon decay signal. The fraction of neutrinos from $\pi^{+}$ decay is 9 (17)% percent in the FHC (RHC) beam. The first population consists of $K^{+}$ at moderate longitudinal momentum, $p_{||}\stackrel{{\scriptstyle>}}{{\sim}}30$ GeV, and with a range of magnitudes of momenta transverse to the proton beam direction, $p_{T}$. The second population consists of $\pi^{+}$ and $K^{+}$ at higher $p_{||}$, $\stackrel{{\scriptstyle>}}{{\sim}}40$ GeV and $p_{T}\stackrel{{\scriptstyle<}}{{\sim}}0.15$ GeV. The second population is dominated by mesons which are underfocused or entirely unfocused by the horns and are common to the FHC and RHC predictions, whereas the first is a unique contribution in the FHC beam since these high-$p_{T}$ $K^{+}$ are defocused in the RHC beam.

Neutrino interactions are simulated using the GENIE neutrino event generator [17] version 2.12.6. Quasi-elastic (1p1h) interactions are simulated using the Llewellyn-Smith formalism [18] with the vector form factors modeled using the BBBA05 model [19]. The axial vector form factor uses the dipole form with an axial mass of $M_{A}=0.99$ GeV/c². Resonance production is simulated using the Rein-Sehgal model [20] with an axial mass of $M_{A}^{RES}=1.12$ GeV/c². Higher invariant mass interactions, including Deep Inelastic Scattering (DIS), are simulated using a leading-order pQCD model with the Bodek-Yang prescription [21] for the modification at low square of the momentum transfer, $Q^{2}$.

A relativistic Fermi gas model [22] with an additional Bodek-Ritchie high momentum tail [23] is used to describe the nuclear environment. The maximum momentum for Fermi motion is assumed to be $k_{F}=0.221$ GeV/c. GENIE models intranuclear rescattering, or final state interactions (FSI), of the produced hadrons using the INTRANUKE-hA package [24].

To better describe MINERvA data, a variety of modifications to the interaction model are made. To better simulate quasielastic events, the cross section is modified as a function of energy and three momentum transfer based on the random phase approximation (RPA) part of the Valencia model [25, 26] appropriate for a Fermi gas [27, 28]. Multi-nucleon scattering (2p2h) is simulated by the same Valencia model [29, 30, 31], but the cross section is increased in specific regions of energy and three momentum transfer based on fits to MINERvA  data [32] in a lower energy beam configuration. Integrated over all phase space, the rate of 2p2h is increased by 50% over the nominal prediction. Based on fits done in Ref. [33], we decrease the non-resonant pion production by 43% and reduce the uncertainty compared to the base GENIE model uncertainties. This modified version of the simulation is referred to later in this paper as MINERvA Tune v1.

The response of the MINERvA detector is simulated using GEANT4 [15] version 4.9.3p6 with the QGSP_BERT physics list. The optical and electronics performance is also simulated. Through-going muons are used to set the absolute energy scale of minimum ionizing energy depositions by requiring the average and RMS of energy deposits match between data and simulation as a function of time. A full description is found in Ref. [13]. Measurements using a charged particle test beam [34] and a scaled-down version of the MINERvA detector set the absolute energy response to charged hadrons. The effects of accidental activity are simulated by overlaying hits in both MINERvA and MINOS from data corresponding to random beam spills appropriate to the time periods in the simulation.

A charged-current $\nu_{\mu}$ event is selected by matching a reconstructed muon track in MINERvA with a momentum and charged analyzed muon track in MINOS. The approximately 252(132) event IMD sample, prior to $E_{\mu}$ and $\theta_{\mu}$ selections, in FHC(RHC) is a small subsample, approximately $0.006\%$ of the inclusive $\nu_{\mu}$ charged-current sample with a reconstructed neutrino interaction point in the tracking fiducial volume. For the high-energy muons in the IMD sample, the selection of $\mu^{-}$ using the direction of the bend in the magnetic MINOS spectrometer is 99% efficient, and in the RHC sample where most muons are $\mu^{+}$, the purity for $\mu^{-}$ selection is over 97%.

Since IMD only produces an energetic forward muon in the final state, the visible energy in the tracker and calorimeters is expected to come only from the muon. By contrast, events from background reactions on nuclei almost always produce some visible recoil, including a recoiling target nucleon. The IMD selection requires visible hadronic energy in tracker and electromagnetic calorimeters be less than $80$ MeV and visible energy within $150$ mm of the neutrino interaction point be less than $10$ MeV. This effectively removes most events with low-energy protons or pions in the final state while retaining all but 4% of the IMD events.

Figure 3 shows the expected and observed signal sample as a function of reconstructed muon energy in the FHC and RHC beams after these selections. The expected signal and background are nearly comparable, and the backgrounds are almost entirely due to charged-current quasielastic scattering, $\nu_{\mu}n_{\rm bound}\to\mu^{-}p$, with a small fraction of events from the multi-nucleon version of this scattering, the 2p2h process described above.

The prediction of the MINERvA flux [12, 6, 7] described in Section II gives a nominal flux prediction, the “central value”, and a series of flux “universes” that describe the uncertainties and covariances in those uncertainties by the Monte Carlo method. The consistency of each flux universe, denoted by $\Phi$ with a universe index, $i$, with the number of IMD events, $N$, is measured by the probability of the measurement given this flux, $P(N|\Phi_{i})$. Since this measurement consists of two weakly correlated measurements in the two beams, the consistency of these measurements with a given flux universe is given by

where $\Delta$ is a vector of the difference between the number of measured and predicted events in the FHC and RHC beams, and $V$ is the covariance matrix of these measurements. Figure 7 shows the predicted number of IMD events in the FHC and RHC beams with the measurement superimposed. It is evident from this that some flux universes are significantly less consistent with the measurement than others. The a priori prediction of the flux can then be modified, according to Bayes’ Theorem, by weighting the flux universes by the probability given in Eq. 4 when forming the central value prediction or the variance of the ensemble of universes. The evident correlation between the FHC and RHC predictions is in part due to the common source of high-energy $\nu_{\mu}$ in the two beams of unfocused low $p_{T}$ parent mesons, as discussed in Sec. II.

The predicted RHC and FHC fluxes and their uncertainties, before and after the IMD constraints are applied, are shown in Fig. 8. As expected, the IMD sample provides a significant constraint for the highest energy neutrinos in the NuMI beams. The flux is modified and constrained at neutrino energies below the threshold of $\approx 11$ GeV because high-energy mesons that decay to make these neutrinos may also decay at larger angles with respect to the beam axis to produce lower energy neutrinos. Therefore, even though we only measure the IMD rate above threshold, we are constraining flux universes that encode the physics that leads to that high-energy part of the neutrino spectrum. The integral flux above $11$ GeV in the FHC beam is predicted to be $2.61\times 10^{-5}\pm 2.28\times 10^{-6}$ $\nu_{\mu}$/POT/m² before the constraint and is evaluated as $2.38\times 10^{-5}\pm 1.50\times 10^{-6}$ $\nu_{\mu}$/POT/m² after.

The MINERvA experiment has successfully isolated a sample of inverse muon decay events, $\nu_{\mu}e^{-}\to\mu^{-}\nu_{e}$, and has used those events to constrain the flux of high-energy neutrinos in its beam. The constraint provides an in situ way to reduce uncertainties from its high-energy flux. Such a method can be applied to any accelerator neutrino beam produced by protons of energies much greater than the $11$ GeV threshold for inverse muon decay, and in particular can be used for a similar purpose in the planned DUNE experiment.