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

Inverse muon decay, $\nu_\mu e^-\to\mu^-\nu_e$, is a reaction whose cross-section can be predicted with very small uncertainties. It has a neutrino energy threshold of $\approx 11$ GeV and can be used to constrain the high-energy part of the flux in the NuMI neutrino beam. This reaction is the dominant source of events which only contain high-energy muons nearly parallel to the direction of the neutrino beam. We have isolated a sample of hundreds of such events in neutrino and anti-neutrino enhanced beams, and have constrained the predicted high-energy flux.


I. INTRODUCTION
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 * Now at Brookhaven National Laboratory † Now at Los Alamos National Laboratory ‡ Now at University of Minnesota 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 arXiv:2107.01059v3 [hep-ex] 23 Nov 2021 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 νe − → ν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), ν µ e − → ν e µ − . The IMD process has a threshold energy of E min = m 2 µ −m 2 e 2me ≈ 11 GeV, and a total cross section given at tree level by [9] where m µ,e are the masses of the muon and electron, G F is the Fermi constant, and the relativistic invariant quantity, s, is the square of the center-of-mass scattering energy. When E ν is measured in the lab frame, s = 2E ν m e + m 2 e . The spectrum of muons emitted for a fixed neutrino energy in the lab frame, E ν , is approximately uniform with limits between E min and E ν , with small corrections to the uniformity and the kinematic limits of order m e /E ν and m e , respectively. Radiative corrections to the process have been calculated, and these decrease the tree level prediction above by several percent, with the largest decreases at the lowest neutrino energies and the kinematic limits [9]. The kinematics of IMD require where y ≡ E µ /E ν and θ µ is the muon angle with respect to the incoming neutrino direction. Practically, this means the muon will be very close in direction to the incoming neutrino. There is a related inverse muon decay process,ν e → e − → µ −ν µ , with identical kinematics and a practically indistinguishable final state. In our experiment, the number ofν e above threshold is at most a few percent of the number of ν µ above threshold, so this contribution is unimportant.
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 ν µ 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 ν µ 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 µ and higher θ µ than the IMD signal.

II. THE MINERVA DETECTOR AND SIMULATION
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 ν µ in the FHC and RHC beams, particularly those originating from π + 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 ± 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×10 20 protons on target (POT) in the FHC configuration and 11.24 × 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 π + decay is 9 (17)% percent in the FHC (RHC) beam. The first population consists of K + at moderate longitudinal momentum, p || > ∼ 30 GeV, and with a range of magnitudes of momenta transverse to the proton beam direction, p T . The second population consists of π + and K + at higher p || , > ∼ 40 GeV and p T < ∼ 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. Quasielastic (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 2 . Resonance production is simulated using the Rein-Sehgal model [20] with an axial  [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.

III. SELECTION OF INVERSE MUON DECAY EVENTS
A charged-current ν µ 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 µ and θ µ selections, in FHC(RHC) is a small subsample, approximately 0.006% of the inclusive ν µ 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 µ − using the direction of the bend in the magnetic MINOS spectrometer is 99% efficient, and in the RHC sample where most muons are µ + , the purity for µ − 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.
In addition, the reconstructed muon must be a µ − with total energy greater than 10 GeV, a threshold below the kinematic threshold of the IMD process due to the ≈ 11% fractional energy resolution in MINOS. Negatively charged muon candidates are also required to have reconstructed energy below 50 GeV, a point beyond which the charge of the muon cannot be reliably measured. From Eq. 2, the kinematics of neutrino scattering from atomic electrons requires that E µ sin 2 θ µ = 2m e 1 − Eµ Eν , where θ µ is the scattering angle with respect to the initial neutrino direction. In a given event, we measure E µ and θ µ , but we do not know E ν a priori, nor do we measure θ µ with sufficient precision to extract E ν from the relationship above. However, we have a minimum E µ for IMD events, and we have a maximum relevant E ν set by our flux which falls steeply with energy as shown in Fig. 1. These facts together imply that

1−
Eµ Eν will typically be a number significantly less than 1, thus allowing us to place a tighter selection on E µ sin 2 θ µ than the maximum of 2m e which is reached in the limit of Eµ Eν → 0. We define E max ν ≡ 35 GeV, and form where the small angle expansion has been used to set sin θ µ ≈ θµ 1radian . The event selection then requires F (E µ , θ µ ) < 2m e . This cut will be increasingly inefficient for neutrino energies above E max ν , but the choice of 35 GeV is predicted to include 98% of the IMD events in the FHC beam and 75% in the RHC beam before accounting for experimental resolutions. The distribution in F (E µ , θ µ ) after all selections and background tuning is shown in Fig. 5. 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 chargedcurrent quasielastic scattering, ν µ n bound → µ − p, with a small fraction of events from the multi-nucleon version of this scattering, the 2p2h process described above.

A. Background Constraints
To constrain the remaining background a sideband sample is measured using events which pass the recoil and vertex energy criteria, but have muon energy between 7 and 9 GeV. The sample composition is almost exactly the same as the backgrounds in the signal selection, but this sample has almost no signal component. Figure 4 shows the sideband sample as a function of F (E µ , θ µ ). As can be seen particularly with the FHC sample, there are two differences between the sideband simulation and data, both likely due to poorly modeled nuclear effects. The first is the overall rate, which will be strongly affected by the probability that outgoing nucleons reinteract to produce neutrons in the final state, which in turn go undetected and allow the events to pass the recoil cuts. The second is that the events at low F (E µ , θ µ ) are suppressed, possibly due to Pauli blocking or nuclear screening. This sideband sample is used to divide the background into events with F (E µ , θ µ ) below and above 2m e , where the former is the one that directly enters into the background subtraction for the signal sample. However, the absolute data and simulation differences between the sideband and the signal with E µ > 10 GeV are also affected by uncertainties in the flux itself. The flux uncertainties are rapidly changing in the sample since the focusing peak is at E ν ∼ 7 GeV. Therefore a second sample of events with E µ > 10 GeV but with 4 ≤ F (E µ , θ µ ) < 10 MeV is added to provide an absolute normalization to the background prediction. The resulting scale factors and their uncertainties are shown in Table I, and the corrected simulated distributions compared to the sideband data are shown in Fig 4. The net effect is to increase the backgrounds compared to the prediction by approximately 15% in the high F (E µ , θ µ ) region, but to suppress the background in the signal region of F (E µ , θ µ ) < 2m e .
After the sideband fit, the scale factors are applied to the selected signal sample. The resulting distribution is shown as a function of F (E µ , θ µ ) and muon energy in Figs. 5 and 6, respectively.

B. Systematic Uncertainties
Systematic uncertainties in this analysis fall under three different categories: flux, detector response, and neutrino interaction model uncertainties. The uncertainties from individual sources are evaluated by reextracting background subtracted samples using modified simulations. The size of each modification is related to the uncertainty in each source. Neutrino interaction model uncertainty in the result is solely due to the background interactions since the signal interaction model is well known.
The flux uncertainty is a typical leading uncertainty in neutrino cross section measurements, but since in this analysis the output is just a count of the number of events, the flux uncertainty enters only through the background constraint which is extrapolated from lower muon and presumably neutrino energy, to higher energy. The resulting small uncertainties from the input flux on the background subtraction are compared with the a priori flux when this result is applied as a flux constraint. A related uncertainty which must be factored into the predicted number of events is the normalization uncertainty of 1.4% from uncertainty in the number of electrons in the target, based on material assays and weight measurements of scintillator planes.
The uncertainty in the detector response to hadrons is evaluated using shifts determined by in situ measurements of a smaller version of the detector in a test beam [34]. Uncertainties in inelastic interaction cross sections for particles in the detector material are independently varied based on data-Monte Carlo differences between GEANT particle cross sections and world data on neutrons [35][36][37][38], pions [39][40][41][42], and protons [43][44][45]. The muon reconstruction uncertainty is dominated by uncertainty in the energy scale, which is constrained by a combination of data and simulation described in Ref. [46] to 1.0%. The uncertainty in the matching efficiency is from imperfect modeling of the efficiency loss from accidental activity in the MINOS near detector when matching muon tracks from MINERvA to MINOS. This last efficiency is also determined by a data-simulation comparison as a function of instantaneous neutrino beam intensity. The interaction model uncertainties are evaluated using the standard GENIE reweighting infrastructure with additional uncertainties from MINERvA Tune v1. The sideband constraint reduces those uncertainties by more than a factor of two.
The final samples have 127(56) selected events in data for the FHC(RHC) configurations. Due to the limited size of the sample, each is only reported as total number of events. The statistical and systematic uncertainties in the background subtracted samples are shown in Tab. II. In the both beams, the prediction is larger than the observed number of events as shown in Fig. 7. Both results are dominated by statistical uncertainty with subleading contributions from the uncertainties in the interaction cross section model and the muon reconstruction. Thē ν e initiated reaction described above is predicted by the cross-section in Ref. [47] to be 0.5%(2%) of the signal rates above, and for convenience was treated as a background in this analysis.

IV. FLUX CONSTRAINTS FROM IMD
The prediction of the MINERvA flux [6,7,12] 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 Φ 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 |Φ 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 ∆ 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 ν µ 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 ≈ 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 × 10 −5 ± 2.28 × 10 −6 ν µ /POT/m 2 before the constraint and is evaluated as 2.38 × 10 −5 ± 1.50 × 10 −6 ν µ /POT/m 2 after.

V. CONCLUSIONS
The MINERvA experiment has successfully isolated a sample of inverse muon decay events, ν µ e − → µ − ν e , and has used those events to constrain the flux of highenergy 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.