Neutrino-electron elastic scattering for flux determination at the DUNE oscillation experiment

We study the feasibility of using neutrino-electron elastic scattering to measure the neutrino flux in the DUNE neutrino oscillation experiment. The neutrino-electron scattering cross section is precisely known, and the kinematics of the reaction allow determination of the incoming neutrino energy by precise measurement of the energy and angle of the recoiling electron. For several possible near detectors, we perform an analysis of their ability to measure neutrino flux in the presence of backgrounds and uncertainties. With realistic assumptions about detector masses, we find that a liquid argon detector, even with limitations due to angular resolution, is able to perform better than less dense detectors with more precise event-by-event neutrino energy measurements. We find that the absolute flux normalization uncertainty can be reduced from ~8% to ~2%, and the uncertainty on the flux shape can be reduced by ~20-30%.


I. INTRODUCTION
The Deep Underground Neutrino Experiment (DUNE) is designed to measure CP violation in neutrino oscillations by making precise measurements of the neutrino flavor oscillations ν µ → ν µ and ν µ → ν e , and their antineutrino analogues, as a function of the neutrino energy, E ν [1]. DUNE uses a wideband neutrino beam peaked at 2.5 GeV, and with 92% of the muon-neutrino flux in the energy range 0.5-5 GeV [2]. The DUNE far detector (FD) will measure neutrino-argon interactions, and infer the neutrino energy from the observed final state particles [3]. In addition to oscillation parameters, these measurements are sensitive to several inputs, each of which has significant, O(10%), uncertainties: the neutrino-argon interaction cross sections, the relationship between the true and inferred E ν , and the neutrino flux. To achieve its physics goals, in particular the measurement of CP violation, the DUNE near detector (ND) must constrain the uncertainties on the predicted event spectra to the level of ∼2-3% [1].
Neutrino cross-section uncertainties are energy dependent, and affect both the rate of interactions and the energy reconstruction. Near detector measurements of neutrino-argon charged-current interactions are extremely valuable, but typically constrain a product of flux and cross sections. The near and far detectors see different neutrino fluxes, primarily because of oscillations, which limits the ability to extrapolate these ND constraints to the FD. The flux as a function of E ν is a priori poorly predicted, primarily because of hadron production uncertainties as described in Section IV. This makes it difficult for the near detector to simultaneously measure the flux of neutrinos and to study the mechanisms by which neutrinos interact.
A helpful way to break this degeneracy is to separately measure the flux of neutrinos as a function of energy at the near detector. This can be done by selecting a sample of events for which the cross section as a function of energy is known. For example, at E ν 1 GeV, the neutrino interaction cross section for events with low energy transfer, ν, is roughly constant with neutrino energy. This way to measure flux, referred to as the "low-ν" technique, has been used to study total and deep-inelastic neutrino cross sections as a function of neutrino energies [4][5][6][7][8][9]. However, the DUNE first and second oscillation maxima occur at ∼2.5 GeV and ∼0.8 GeV, respectively. At these energies, it is difficult to select a sample with ν/E ν sufficiently small, and the low-ν technique breaks down. In this energy range, neutrino-electron elastic scattering is the only process with a known cross section.
In this paper, we demonstrate the feasibility of measuring the neutrino flux as a function of energy in the DUNE near detector using neutrino-electron elastic scattering. The LBNF beamline simulation is used to produce a flux of neutrinos at the near detector location, including the effect of the beam dispersion, as described in Section IV. We use Geant4 simulations to study the expected electron angular resolution in the liquid argon time projection chamber (LAr TPC) of the DUNE ND in Section V. Similar studies are performed for a highresolution gaseous detector, and a plastic scintillator detector. A detector with perfect electron reconstruction and background rejection is also considered as a limiting case. Section VI describes the details of the flux fits, and presents the results. In Section VII, we discuss the potential for using inverse muon decay to further constrain the high energy ν µ flux. Finally, in Section VIII, we present our conclusions.

II. NEUTRINO ELECTRON SCATTERING
Neutrino-electron elastic scattering, νe − → νe − , is precisely predicted by the electroweak theory because it is a 2 → 2 process that involves only weak interactions of fundamental leptons. In the limits that the neutrino energy E ν is much greater than the electron mass m e and far below the energies required for resonant W boson production, E ν M 2 W 2me , the νe − → νe − cross section for neutrinos or antineutrinos is given at tree level by Here, G F is the Fermi weak coupling constant, s is the Mandelstam invariant representing the square of the total energy in the center-of-mass frame, and y ≡ T e /E ν where T e is the electron kinetic energy. The couplings C LL and C LR are different for neutrinos and antineutrinos and depend on flavor. For ν µ and ν τ , C LL = − 1 2 + sin 2 θ W and C LR = sin 2 θ W , where θ W is the Weinberg angle, and in the corresponding antineutrino couplings, the values for C LL and C LR are interchanged. For ν e (ν e ), the value of one of the couplings, C LL (C LR ), is 1 2 + sin 2 θ W because of interfering contributions from neutral-current interaction that is present for all flavors and from a chargedcurrent interaction that is present only for electron neutrinos. Electroweak radiative corrections to the process are few-percent corrections and are discussed in detail in Appendix A.
The theoretical uncertainty of the neutrino-electron elastic scattering cross section from uncertainties in the parameters and radiative corrections is small [10]. Recent work [11] has shown that the limiting uncertainty comes from hadronic loops in radiative corrections which results in a few permille uncertainty. Therefore a measurement of the reaction can be used to measure neutrino flux at this precision. At the ∼ O(1) GeV neutrino energies of DUNE, this cross section is approximately 10 −4 of the total charged-current ν µ cross section; therefore the number of events is small and backgrounds may be substantial. However, for realistic near detector sizes, the event sample is expected to be sufficiently large in the DUNE beam to allow for statistical precision on a neutrino-electron elastic scattering sample to be O(1%) [12].
The angle of the final state electron with respect to the neutrino, θ e , is where E e is the energy of the final state electron. Therefore at neutrino energies ∼1 GeV, such as for DUNE, where m e E ν , the final state electron is very forward. A measurement of the angle and electron energy determines y, and thus also the neutrino energy.
Another neutrino-electron scattering process with a well-known cross section is inverse muon decay (IMD), ν µ e − → ν e µ − . This process has a threshold energy of ≈ 11 GeV, and a total cross section given at tree level by [13] 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. This cross section increases with energy as the DUNE flux is falling, and the event rate is expected to peak at ∼18 GeV. IMD could provide a constraint on the high energy tail of the ν µ flux; however, such a constraint would have little impact on the DUNE neutrino oscillation analyses. This process is discussed further in Section VII, and in less detail than νe − → νe − in this manuscript.

III. MINERVA'S NEUTRINO-ELECTRON SCATTERING FLUX MEASUREMENT
The MINERvA experiment is the only accelerator experiment to date that has successfully used this technique [14,15] to significantly reduce its uncertainty on a predicted neutrino flux. MINERvA reconstructed these events in a segmented scintillator detector with neutrinos at energies similar to DUNE's. The first analysis with the low energy NuMI beam [14] observed 127 total events including a predicted background of 30 ± 4 events; a second, recent analysis with the medium energy NuMI beam [15] found 1021 events with a predicted background of 212 ± 13.5 events. The background composition of the two analyses was different because of the event selection and beam energies. In the medium (low) energy analysis, the background was approximately 28% (55%) ν e charged-current interactions, primarily quasielastic like events ν e n bound → e − p, 54% (30%) neutral current interactions, primarily with a π 0 in the final state, and 18% (15%) ν µ charged current events, also primarily with a π 0 in the final state and a very low energy final state µ − . In both analyses, backgrounds in the segmented scintillator were reduced by requiring an electron energy of 800 MeV or greater, which is not a desirable selection for a DUNE near detector because of the physics interest in the low energy neutrino flux. Because of the angular resolution in the MINERvA segmented scintillator, with a granularity of ∼2 cm, MINERvA did not attempt to use angular information to reconstruct the incoming neutrino energy. The systematic uncertainty on the observed rate in the MINERvA medium (low) energy measurement was 1.8% (5%), and was mostly due to uncertainties in the background reactions. The uncertainty on background reactions, particularly the low Q 2 behavior of the ν e quasielastic-like background events, is significantly lower in the medium energy analysis than in the low energy analysis, largely due to better knowledge of the low Q 2 behavior of neutrino reactions due to MINERvA data itself [15,16]. In the medium energy analysis, the electron reconstruction efficiency and electromagnetic energy scale of the detector were also noted contributors of systematic uncertainty, but were not dominant sources. Both analyses had a 1-2% uncertainty in the application of the event rate to the neutrino flux prediction from the fiducial mass of the detector. As a result of these analysis, the fractional uncertainty on MINERvA's low energy flux between 2 and 10 GeV was reduced from 8.7% to 6.0%, and the uncertainty at the focusing peak of the medium energy beam was reduced from 7.5% to 3.9%.

IV. LBNF BEAM
The planned Deep Underground Neutrino Experiment (DUNE) [1,3] will operate in the Long Baseline Neutrino Facility (LBNF) [2] beamline at Fermilab. LBNF is designed to operate at an initial beam power of 1.2 MW, with a design capacity of 2.4 MW, more than three times the maximum intensity of the NuMI beamline (700 kW) [17]. At 1.2 MW intensity, corresponding to 1.1×10 21 protons on target per year, and with a detector located 574m from the neutrino source, ∼120 ν-e − events are expected per year per ton of argon. Hydrocarbon detectors have a higher ratio of electrons to nucleons, and therefore have higher event rates per unit mass. The expected rates per year per ton are ∼144 for CH and ∼152 for CH 2 .
To create the LBNF neutrino beam, protons from the Main Injector strike a fixed target, producing pions and kaons, which are focused by a system of magnetic horns into a decay volume, where the mesons decay primarily into muons and muon-flavor neutrinos. The specific parameters of LBNF are optimized to maximize the sensitivity to CP violation in DUNE. The target is a long, thin rod, 2 m in length and 10 mm in diameter. Three horns are used in the optimized configuration of LBNF, compared to only two horns in NuMI. The main advantage of the third horn is a reduction in the high-energy tail of the flux. The target protrudes slightly into the first horn, while the second and third horns are located 3 m and 17 m downstream of the target, respectively. The decay pipe volume is cylindrical, 200 m along the axis, and with a radius of 2 m. It begins 20 m downstream of the target, and is angled downward at 6 degrees (101 mrad.) with respect to the horizontal, such that it points toward the on-axis far detector.
The horns are designed to focus positive mesons in forward horn current (FHC) mode, leading to a primarily ν µ beam. In reverse horn current (RHC) mode, the beam consists primarily ofν µ . The wrong-sign contamination is higher in RHC because the proton-carbon interactions produce more π + than π − , but the contamination is primarily in the flux tail. Electron-flavor neutrinos make up 1% of the total flux, and arise primarily from the decay chain π → µ → e. At energies above 10 GeV, neutrinos from kaon decays dominate, and the ν e contamination is larger from K ± → π 0 e ± ν e and K 0 L → π ± e ∓ ν e . The neutrino flux is peaked at 2.5 GeV, the oscillation maximum for a baseline of 1300km.
The beamline is simulated with g4lbnf, a Geant4-based model. Proton interactions in the target, as well as the subsequent interactions of hadrons in the target and focusing system, are simulated using the "QGSP_BERT" physics list, which combines the quark-gluon string with precompound (QGSP) model and a Bertini cascade at higher energies. This analysis is based on g4lbnf version v3r5p4, which is based on an 120 GeV proton beam. The energy spectra for all four neutrino species in both horn polarities used in this analysis are shown in Figure 1.
Neutrinos can have non-zero angles with respect to the beam axis due to imperfect focusing, the finite width of the decay pipe, and the finite size of the detector. The mean neutrino angle at the LBNF near detector facility, 574m from the target along the beam axis, is approximately 1.5 mrad. The decay pipe geometry gives a maximum angle of 5.6 mrad. to the center of the near detector, corresponding to a decay at the edge of the decay pipe and nearest to the detector hall.
Larger angles are possible for neutrinos originating from decays outside the decay pipe region (88% of neutrinos intersecting the near detector originate in the decay pipe, and an additional 8% in the target hall). Muon decays can produce neutrinos with much shorter baselines. The area-normalized angular distribution for FHC is shown in Figure 2 for the four neutrino species. The antineutrino distributions are more sharply peaked at zero angle because of very forward pions passing down the center of the horns, which are not defocused. These pi- ons are generally higher in energy and give very forward decays. The larger high-angle tails for antineutrino and ν e are due to muon decays. The reconstruction of the neutrino angle is critical in an analysis of νe − → νe − . The angle due to the finite detector width can be corrected by taking the neutrino angle to be the line connecting the mean decay position to the reconstructed interaction vertex. In LBNF, the mean decay position is approximately the center of the decay pipe due to approximate cancellation of the exponentially decreasing pion flux and quadratically increasing detector solid angle. The angular distribution due to the decay pipe width cannot be corrected, and effectively smears the distribution of θ e , the electron angle with respect to the neutrino.
The largest uncertainties on the flux prediction as a function of neutrino energy are due to hadron production on the target, and tolerances of the focusing system. The method for evaluating these uncertainties is similar to that of MINERvA [18]. There are strong positive correlations between the flux predictions at the near and far detectors, between forward and reversed horn current, and between muon and electron (anti)neutrinos. Because of these correlations, a near detector flux constraint also constrains the far detector flux. In principle, an FHC (RHC) constraint would also constraint the RHC (FHC) flux, although in this analysis we investigate constraints in both modes.
Uncertainties on the distribution of the incoming neutrino angle could bias an extraction of the energy spectrum from νe − → νe − . The effect of hadron production uncertainties on the angle are found to be negligibly small. Varying the horn focusing parameters gives sub-percent changes to the angular distribution. One exception is the uncertainty on the width of the decay pipe, which determines the endpoint of the distribution in Fig. 2 by specifying the maximum possible angle from the decay pipe to the near detector. This uncertainty produces large changes to the flux at very high angle, but affects less than 1% of the total flux.

V. NEAR DETECTOR TECHNOLOGIES STUDIED
For this study, we evaluate the expected performance of the DUNE near detector for νe − → νe − . The primary design considered is a liquid argon time projection chamber (LAr TPC), roughly based on the ArgonCube concept [19]. We also consider a high-resolution tracking detector (HRT), which has significantly improved energy and angular resolution compared to the LAr TPC, but a lower mass because it would require a gaseous rather than liquid argon target. This is meant to represent a generic tracking detector with a large number of high resolution spatial measurements per unit of dE/dx, but the performance parameters are roughly what could be achieved with a straw tube tracker like the one described in the reference design of the DUNE conceptual design report [3], or with a gaseous argon TPC with good angular resolution. The purpose of including it in the study is to determine whether a lower mass detector with superior energy and angular resolution can provide a stronger constraint. In addition, we consider a solid plastic scintillator (labelled "CH") detector, which is essentially what could be achieved by putting the MINERvA detector [20] in the DUNE beam. In this analysis, we assume a 5 year exposure on each detector, with each horn current (the FHC and RHC analyses are performed separately). The LAr fiducial mass is assumed to be 30 tons, which can easily be achieved with a total active volume of 4 × 3 × 5 m. Alternatively, this can represent a fiducial mass of 60 tons in a scenario where only 50% of the total exposure is taken on-axis. The HRT and plastic scintillating detector are assumed to have a 5 ton fiducial mass. For comparison, the HRT and CH analyses are also repeated with a 30 ton mass equal to the LAr.
There are three detector parameters that impact performance in the ν-e − channel: the electromagnetic energy resolution, the angular resolution for forward electrons, and the threshold for rejecting events with other final-state charged particles, such as low-energy protons. Angular resolution is the most important metric for νe − because the signal kinematic limit and the neutrino energy reconstruction both depend on E e θ 2 e . In this section, the procedure used to determine the expected angular resolution as a function of electron energy is described in detail for the liquid argon ND concept. The procedure is qualitatively similar for the HRT and CH detectors.
The liquid argon concept for the DUNE ND is based on ArgonCube. The detector consists of an array of optically segmented TPC modules in a common cryostat. The module size has a cross section of 1× 1 m 2 , with a central cathode dividing the TPC into two drift regions with a maximum drift distance of 50cm. Charge is read out by an array of pads, instead of the projective wires in the DUNE far detector design. We consider a pad size of 3 × 3mm, similar to what was used in initial demonstrations of the pad readout technology [19], for a total of ∼10 5 channels per m 2 . With a maximum drift length of 50 cm, transverse diffusion is estimated to be 0.8 mm, based on 13 cm 2 /s at 1 kV/cm [21][22][23].
The angular resolution is determined by the position resolution of the detector, and by multiple scattering of electrons in liquid argon. Forward-going electrons from ν-e − elastic scattering are nearly parallel to the readout plane, and will intersect individual rectangular pads at nearly right angles. The 2D electron angle in the plane perpendicular to the drift direction depends on the pad coordinate; for 3mm pitch the position resolution is 3mm / √ 12 = 0.87mm. A potential aliasing effect exists for tracks nearly parallel to a row of pads, which can be mitigated by staggering the pads in successive rows. Improved position and thus angular resolution can be achieved with a triangular pad design. Because diffusion is small due to the short drift, collected charge for a forward-going particle will be shared among two adjacent triangles, and the relative charge collected on each triangle is proportional to the position of the electron. This feature is used by MINERvA to achieve 3mm position resolution with 1.7cm scintillator strip pitch. The 2D angle in the drift plane is determined by timing. The neutrino interaction time is determined by the detection of scintillation photons. The position resolution in the timing direction is expected to be significantly better than the pad size.
The radiation length in liquid argon is 14cm. For a 3mm pad pitch, this corresponds to N = 47 position measurements per radiation length, X 0 . The resolution due to measurement, σ meas , and multiple coulomb scattering, σ MS , can be calculated as where L is the track length. The measurement resolution decreases with track length, while the multiple scattering term increases. The optimal track length to fit is approximately one radiation length. Other hard scattering processes also contribute. To quantify this effect, electrons are simulated in liquid argon with Geant4 10.2. The electron position is determined at 3mm intervals and smeared by a Gaussian function with a sigma of 1 mm. An uncertainty is placed on each position measurement based on the average multiple scattering according to Equation 4. The resulting points are fit to a straight line. The reconstructed angle is then compared to the true angle to determine the resolution. This procedure is repeated for electrons of varying momentum. Figure 3 shows the residual of the twodimensional angle, ∆θ e,x , for electrons at 1 GeV and 5 GeV. The energy dependence is parameterized as a function of electron energy according to a double Gaussian. The widths of each Gaussian, and the relative normalization, are determined from fits to these distributions. The dependence on electron energy is fit assuming constant and 1/E e terms. The central multiple scattering term plateaus at 4 mrad., which is essentially the measurement limit. The dependence on electron energy is shown in Figure 4.
The expected angular resolution functions for the HRT and plastic scintillator detectors are determined similarly. For the HRT, a position measurement is made in each two-dimensional projection every 8 cm, with an assumed transverse resolution of 200 µm. The angular resolution function is determined by the same procedure as described above, but in a volume with density 0.1 g/cm 3 , which is approximately the average density for a strawtube or high-pressure gaseous detector. An electron density equivalent to that of CH 2 is assumed. The resulting angular resolution in the HRT is significantly better than in LAr, reaching ∼1.5 mrad. at high energy. The scintillator detector is based on MINERvA [20]. MINERvA tracker planes are spaced by 2 cm, with at least 4 cm between two planes of the same twodimensional orientation. The assumed position resolution per plane is 3 mm. The angular resolution achieved with these assumptions is somewhat worse than in LAr, ∼8 mrad. at high energy, and consistent with the angular resolution of MINERvA.
For all detector types considered above, the reconstructed electron energy is treated in the same way. It is parameterized by a Gaussian centered on the true electron energy with a width of 5%, and a low-side tail, which affects 10% of events and has a fairly arbitrary 6 form, which serves to smear events in this tail between E true e /2-E reco e without a step function, and includes a generic misreconstruction effect in the analysis. Figure 5 shows an illustrative example, of the E reco e distribution for an electron with E true e = 1 GeV. Systematic studies which vary the size of the low-side tail, and shift the peak of the distribution (to mimic an energy scale bias) are considered in Section VI E.
Finally, for the reference "perfect" detector, the electron energy and angular reconstruction are assumed to be perfect, and perfect background rejection is assumed. The intention with the perfect detector is to show the inherent limitations to the technique due to the relatively low ν-e − event rate and the divergence of the neutrino beam at the near detector location. For the perfect detector, CH 2 was used as the target as it has the highest electron density.

VI. STUDY FRAMEWORK
The aim of this study is to test how well the flux normalization and shape uncertainties can be constrained from the reconstructed electron energy and angle with respect to the nominal beam axis (E e , θ e ) for different potential DUNE near detector designs. We simulate ν-e − scattering events with the GENIE event generator [24,25], which uses the tree level cross section given in Section II and include radiative corrections to the cross section. The full DUNE three-horn optimized flux is used, including the beam divergence as described in Section IV. Detector resolutions, as described in Section V, are used to smear the GENIE prediction as a function of E e , θ e . Additionally, backgrounds from ν e -40 Ar interactions which produce electrons or single photons from π 0 decays are included in the study, as described in Section VI A. The fitting framework is described in Section VI B, the main results are shown in Section VI C, and bias tests to check the robustness of the result to different input assumptions are described in Sections VI D and VI E.

A. Selection and backgrounds
The selection for this analysis is very simple. The signal processes will produce a single very forward-going electron, and no other particles at the vertex. We consider two types of backgrounds: ν e -40 Ar interactions which produce a forward-going electron; and events in which a single photon from a neutrino-nucleus produced π 0 is reconstructed. These are referred to as the ν e and γ backgrounds. We apply a cut on the extra energy deposited at the vertex of E extra ≤ 20 MeV (≤ 30 MeV for the CH detector) for ν e -40 Ar interactions, and a cut of E γ ≤ 50 MeV on the second photon energy for the π 0 production background. Additionally, the π 0 background is suppressed by a factor of 0.1 to account for the γ/e ± separation capabilities of the detectors considered. The extra energy of 20 MeV corresponds to a proton with a range of 5 mm in liquid argon, which will deposit energy on more than one pad. Under these assumptions, the ν e backgrounds are always significantly larger than the π 0 backgrounds. As noted in Section V, for the "perfect" detector options, perfect background rejection is assumed.
Although it was shown in Ref. [14] that a cut on E e θ 2 e provides good separation between signal and background events, we do not make such a cut in this analysis. Instead, we only consider events with a reconstructed electron angle θ reco e ≤ 60 mrad., and perform the fit in E reco e θ reco e space, in which signal and background are reasonably well separated. Because of this separation, and the fit method used (described below), an E e θ 2 e cut is not necessary.

B. Fitting framework
In this analysis, the simulated data, including backgrounds and all resolution effects described in Section V, are binned into E reco e -θ reco e bins, and scaled to the expected event rate given the detector mass and exposure time relevant for each ND options described in Section V. Bins are 4 mrad. wide in θ e , with 15 bins in the range 0-60 mrad. In E e , there are 45 bins in the range 0-60 GeV, where the bin edges are defined such that the central value ±5% lies inside the bin, with a minimum bin width of 0.2 GeV, motivated by the expected E e resolution of ∼5% in our simple model. Although the binning is somewhat arbitrary, changes to the binning (2x finer binning) had no significant effect on the fits described below.
We use a simple template fitting approach, to fit the simulated data with Monte Carlo, as we would for real data. Each template is binned in the same reconstructed θ e -E e bins as the simulated data, and is integrated over a true E ν range. By varying the normalizations of these templates to the fit to the simulated data, a constraint on true E ν can be extracted. The E ν binning is chosen by merging DUNE flux bins such that each merged bin contains a minimum of 500 events, to ensure that the template normalization parameters can be approximated with a Gaussian. The binning for each of the ND scenarios for FHC and RHC is summarized in Table I. Note that each template gives a 2D reconstructed E reco e -θ reco e distribution, which has been integrated over the nominal flux distribution between the E min ν and E max ν boundaries of the template, and over all neutrino flavors (as the data cannot distinguish them), again using the nominal relative fractions of each flavor given by the nominal neutrino flux. These two necessary assumptions are a potential source of bias in the analysis, which will be discussed in detail later. By using an adaptive binning based on the expected statistics, we balance the impact of this bias with the statistical error. Two additional templates are included for the ν e and γ backgrounds, the normalizations for which are also unconstrained parameters in the fit. Example templates are shown in Figure 6 for the 5t CH detector in FHC.
In the fit, we use the "L-BFGS-B" algorithm [26], from the SciPy optimize package [27] to minimize the Poisson-Likelihood: where µ i ( x) is the MC prediction, which is a function of the template normalizations, x, and n i is the number of data events in the ith bin. We exclude bins with E e < 0.5 GeV from the fit to account for a detector threshold.
Modifying the value chosen for the detector threshold had a minimal effect on the fit because there are relatively few events in the very forward, low E e region (which corresponds to very low E ν ), and the lowest E ν template extends well past the threshold in all cases. An example fit, using the nominal LAr design, where a "statistical throw" has been performed on the simulated data, is shown in Figure 7. In each bin, the number of events is drawn randomly from a Poisson distribution; this acts as a very basic sanity check for the fitter. The fitted E reco e -θ reco e distribution approximates the simulated data well. In Figure 8, we show the output correlation matrix from this fit. The "checkerboard" pattern is striking, neighboring bins are strongly anticorrelated, which indicates that the neighboring templates have a very similar effect on the fitted distribution. This is not a problem; indeed, using such fine E ν binning maximizes the flux constraining power of the fit, and minimizes the potential for bias. However, it does mean that the postfit distributions of individual fits are difficult to interpret byeye. Figure 8 also shows that the correlations are small between the signal template parameters, E i , and the two background templates, labeled γ and ν e . This indicates that the fit is able to distinguish signal and background very well, which is unsurprising given the different regions of E reco e -θ reco e space they occupy, as shown in Figure 6. The independence of the signal and background templates will be checked more quantitatively with bias studies in Sections VI D and VI E.
Although the postfit distribution is difficult to interpret by-eye, it can be used to constrain the flux by assigning weights to possible fluxes. The probability of a possible flux being consistent with the measured ν-e −  FIG. 6: Example fit templates for neutrino-electron elastic scattering in six bins of true neutrino energy, and for two background categories, for the 5t CH detector in FHC mode. Each template shows the expected event spectrum as a function of electron energy and angle, for neutrinos in a given energy range.  data can be calculated using [14] P where Σ N is the data covariance (see the correlation matrix in Figure 8), |Σ N | is its determinant, N is the postfit E ν template normalizations, M is the model rebinned to match the template binning, and κ is the number of E ν templates. The probability calculated in Equation 6 can be used to constrain the flux covariance matrix provided by the beam group, ξ ij , to show the impact of the ν-e − constraint. The postfit covariance matrix Ξ ij can be calculated A comparison of the pre-and post-fit covariance matrices can be used to investigate how well the ν-e − sample can constrain the flux. The vector M i give the central values of the post-fit, and deviations of these from the true value are a useful measure of bias in this procedure. Note that because the LBNF beam is a mix of different flavors (see Figure 1), all with different ν-e − cross sections, it is not possible to constrain the flux completely independently of an input flux model, as some assumptions have to be made about the relative contributions from each flavor, if not their spectra. The technique described in this work could be used to constrain any flux model, but we choose to show the additional constraint which can be applied to the Geant4-based DUNE flux simulation, as it is the most sophisticated set of assumptions about the relative contributions from each flavor that we have available.

C. Results
Using the output data covariance from each fit (see the LAr example correlation matrices in Figure 8), and Equations 6 and 7, the constraint on the DUNE flux prediction for each beam and detector configuration can be calculated. In each case, we constrain the full flux covariance, as shown for the nominal LAr configuration for both FHC and RHC in Figure 9, along with the prefit covariance matrices provided by the beam group for comparison. It is clear that the uncertainties are much smaller for the ν µ and ν e (ν µ andν e ) in FHC (RHC) after the fit. The correlations between flavors has also been decreased significantly, although anticorrelations are introduced because the ν-e − sample cannot distinguish flavors, so decreasing the contribution from one flavor can be increased by increasing the contribution from another. Note that the relationship between flavors is already limited by the input beam covariance matrix. The relationship between flavors is more complicated for the RHC case, with stronger correlations obvious in the postfit covariance shown in Figure 9. Although this is expected by the larger wrong-sign contamination in RHC relative to FHC, the fitting technique described here seems to work for both.
Although interesting, the covariance matrices shown in Figure 9 are difficult to interpret by eye, so for the remainder of this work, only the diagonal elements of the covariance will be considered, although we note that the full covariance is calculated each time when producing these plots. The diagonal elements of the covariance matrices are shown for all detector configurations, for both FHC and RHC in Figure 10, with the diagonal elements of the pre-fit flux covariances for comparison. It can be seen that for both FHC and RHC modes, the uncertainty on the dominant flavor (ν µ andν µ respectively) in the flux peak is reduced from ∼8% to ∼1-3%, depending on the detector configuration used.
In order to make it easier to compare the different detectors considered, a ratio is taken with respect to the diagonals of the nominal beam covariance matrix, in Figure 11, and only the dominant ν µ (ν µ ) flavor contributions are shown for FHC (RHC) as they are most interesting. It is clear from Figure 11 that the flux constraining power of the 5t detectors is significantly weaker than for the 30t LAr detector, indicating that the improved reconstruction performance of the CH or HRT does not add significant strength to the flux constraining power of the analysis, except perhaps in the lowest energy bins well below the flux peak. This is also true for the perfect 5t detector, so is not simply a consequence of our assumptions about the HRT/CH performance, it seems that statistics are paramount for this analysis. As expected, the HRT does better than the CH detector, although this may be partially due to the higher electron density in CH 2 than CH. Figure 12 is the same as Figure 11, but with 30t versions of all detector technologies. As expected, with equal masses, the LAr detector performs the least well, which is due to the lower electron density, and worse resolution than the other detector options, but the improvements to the flux constraint seen for the other detector options are fairly small compared with LAr, which supports the conclusion that the most important factor is the statistics gained with larger masses. That said, at a certain point higher statistics would not help, as the intrinsic divergence of the beam would become the limiting factor for the analysis. We have not determined at which LAr or perfect detector mass the limit of statistical improvement occurs in the present study.
By making a shape-only version of the pre-and post-fit flux covariance matrices, by normalizing each flux throw such that the integral is the same, and forming the postfit flux covariance with Equation 7, and then taking a ratio as in Figure 11, the ability for the ν-e − sample to improve the flux shape can be investigated. Such a plot is shown in Figure 13, from which it is clear that there is only a marginal improvement in the understanding of the flux shape relative to the input beam covariance matrix. The shape-only uncertainties are still ∼70% (∼80%) of the nominal shape uncertainties for FHC (RHC). It seems that the power to constrain the flux normalization is largely responsible for the improvements seen in Figures 10 and 11.
Improved detector resolutions would be expected to have a larger impact on the shape constraint than on the total flux normalization, and indeed, we can see from Figure 13 that the 5t perfect detector is nearly comparable to the performance of the 30t LAr detector. In the equal mass case, shown in Figure 14, we see that the 30t better resolution detectors do substantially better than the LAr, but that even for these large detectors, the shape uncertainties are still 50% of their nominal.
In addition to the postfit covariance, another inter-esting quantity is the weighted average flux values, M i , obtained when calculating the covariance matrix in Equation 7. These are shown for all detector configurations considered in both FHC and RHC modes in Figure 15. Large deviations from the pre-fit flux would indicated a bias in some E ν region; no such deviations are observed. More sophisticated studies of possible bias in the procedure are performed in Section VI D.

D. Bias tests
Some bias in the fit results is expected for two reasons. Firstly, the E ν binning used in the fit is coarser than the binning of the flux covariance matrix, at least for some regions of E ν , so variations within those bins cannot be correctly handled by the fit and will introduce small biases. However, the E ν binning is limited by the statistics of the ν-e − scattering sample, so biases due to the coarse binning should be small relative to the statistical uncertainty in the fit. Secondly, the fit cannot distinguish between flavors, and all flavors contribute to the ν-e − sample, but with different cross sections (see Section II). The fit therefore implicitly relies on the expected relationship between flavors, and throws of the flux covariance matrix which change this relationship cannot be dealt with perfectly by the fit. Indeed, one of the weaknesses inherit in using a ν-e − sample as a flux constraint for an accelerator neutrino beamline is that assumptions have to be made about the relative contribution of each flavor.
Although the weighted average flux values obtained when calculating the covariance matrix in Equation 7, and shown for the nominal detector configurations in Figure 15 should give an indication of any bias, and the results suggest that any bias is relatively small, the size of the bias can be calculated more quantitatively by modifying the "data" in the fits, and fitting with the nominal Monte Carlo simulation. If there is an implicit bias towards the input MC, then the fit will not reproduce these modified data distributions. The level of agreement, including uncertainties, between the best fit flux distribution and the input "true" distribution can be assessed with the test-statistic where, the indices i and j are over the true E ν bins used in the fit, and the matrix M ij is the post-fit covariance matrix between fit parameters (as in Figure 8), with the background parameters removed.
To provide a meaningful point of reference, the bias in the fitting technique described in this work was assessed using throws of the input flux covariance matrix. Throws of that matrix, produced through Cholesky de-     composition, describe the expected variation in the neutrino flux given all prior known information about the beam. Figure 16 shows the distribution of χ 2 -values, calculated with Equation 8, obtained for 10,000 independent throws of the input covariance matrix. It is clear from Figure 16 that the biases seen for both FHC and RHC are indeed small for all detector configurations tested. Biases are larger in RHC than FHC, which is probably due the relative beam purities. The biases are larger (most noticeable in RHC) for the 5t detectors than the 30t LAr detector because of the small number of templates used in the fits (see Table I), which was confirmed by checking that the bias disappears for the 30t versions of those detectors. Improvements might be seen by modifying the fit to not expect that fit parameters should be Gaussian, which would relax the 500 expected events per bin requirement which set the template binning. However, reducing the required number of events per bin to 300 did not significantly change the result, so any improve-ment to alleviate that small bias is outside the scope of this work. This bias can be understood as being due to the width of the templates, and the inability of the fit to deal with changes to the assumed flux distribution within each template.
In order to test how sensitive the fit is to changes in relative contribution from each flavor, the post-fit covariance matrix and bias tests were reproduced using a modified version of the input beam covariance matrix, where the covariances between flavors was removed. The preand post-fit covariance matrices, are shown in Figure 17, which can be compared with Figure 9, which included the flavor correlations. The distributions of χ 2 -values for each detector type, calculated with Equation 8, obtained for 10,000 independent throws of the modified covariances matrices from Figure 9 are shown in Figure 18. The biases have uniformly increased slightly with respect to those shown in Figure 16, indicating that the fit relies on the expected flavor composition of the beam as ex- pected, although the bias in the fitted distributions are not significant if those assumptions no longer hold true. Although this test shows a small bias, there is still an implicit reliance on the shape of each flavor contribution to the flux. It should be noted that if the correlation between flavors and the bin-to-bin correlation between each flavor were removed, the fit would not perform well at all because it would have no way to break the ambiguity between the different flavor contributions to the rate.
The bias studies shown in Figure 16 show that the fit converges on the true input flux in an unbiased way for the variations expected given our prior understanding of the beam. However, as well as reducing the flux uncertainty, we would hope that the ν-e − sample could provide an independent check that the flux is correctly modeled, and correctly fit, and identify, distortions to the flux distribution in data which are not covered by the uncertainties in the input flux covariance matrix. To that end, a "crazy" flux distribution was produced independently for each mode by increasing the target density by 30%, well outside its tolerance. Figure 19 shows how well the fit performs when trying to reproduce the "crazy" fake data sets in both FHC and RHC modes. The weighted average flux values obtained by all of the detector configurations tested roughly reproduces the crazy flux reweighting as a function of E ν (indicated by the black line in Figure 19). Although the agreement is not perfect, the post-fit distributions follow the distorted flux well for all detector options in both FHC and RHC mode, indicating that νe − samples would be able to diagnose a large flux bias in both modes, which is reassuring. The ability to fit out the flux distortion worsens at high energies, which is likely to be due to the sparse E ν templates used in the fit (see Table I), but is good across the bulk of the DUNE flux distribution. This constitutes another important result and highlights the necessity of making ν-e − scattering measurements at DUNE. Even though the shape of the flux uncertainty cannot be significantly reduced from the uncertainty obtained from the beam simulation, it can diagnose unmodeled issues with flux prediction in situ.

E. Systematic uncertainties
In the eventual application of this analysis to DUNE, various systematic uncertainties will need to be included or marginalized over in the fit. These will include detector uncertainties, signal and background cross section uncertainties. In order to quantify the effect that such uncertainties have on the analysis, and the size of the uncertainties on the post-fit flux distribution, a series of fake data studies have been performed. In these studies, the Monte Carlo remains unchanged as described above, but instead of using the same Monte Carlo simulation as the "fake data" in the fit, the Monte Carlo is modified to form an alternative fake data set, and the fit is repeated to quantify the effect that this unknown effect has on the output. We investigate a number of systematic sources, described in this section, using the nominal LAr detector configuration.
Here we investigate two sources of systematic error in the energy reconstruction. Firstly, we investigate whether mismodeling the fraction of the reconstructed electron energy in the low-energy tail significantly biases the results. In this analysis, the nominal fraction of events in the low energy tail is 10%, and we have assigned a very conservative ±5% systematic uncertainty. In Figure 20, the error band due to that systematic uncertainty is show on the weighted average flux values (the values of M i produced by Equation 7), relative to the nominal case (a line at y = 1), for the nominal LAr detector in both FHC and RHC. The second energy reconstruction systematic is a 2% uncertainty in the reconstructed electron energy. Again, in Figure 21, the error band due to that systematic uncertainty is show on the weighted average flux values, for the nominal LAr detector in both FHC and RHC. Both, conservative, energy reconstruction uncertainties considered introduce percent-level changes in the average flux value across most of the energy range, and subpercent-level changes around the flux peak, which are small relative to the size of the postfit flux uncertainty (see Figure 10 for comparison) in both FHC and RHC modes. We note also that the variations shown here did not significantly change the average χ 2 value obtained in the 10,000 flux throws relative to the nominal case shown in Figure 16, indicating that although there are shape variations to the weighed average flux distribution in Figures 21 and 20, these are within the expected postfit shape uncertainty in the flux. For larger shifts to the energy reconstruction systematics, larger biases are seen, as would be expected.
Although the beam pointing uncertainty is included in the input flux covariance matrix, it is interesting to ask what the effect on the analysis would be if the beam direction were mismodeled by some constant amount. Fake data studies were performed where a bias of 2 mrad. was introduced in the x-axis of the beam direction, by shifting the reconstructed electron angle w.r.t the nominal beam direction in the MC used in the fit. Figure 22 shows the weighted average flux values for the nominal LAr detector in both FHC and RHC modes. Here, no error band is produced as a bias of ±2 mrad. in any direction would have the same effect as the response is symmetric around the beam axis (at least for the rather simple detector reconstruction considered here). The effect from the beam pointing bias is small for both beam modes. As for the energy reconstruction systematics considered, the beam pointing uncertainty did not increase the average χ 2 for the 10,000 flux throws considered, relative to the nominal case shown in Figure 16, as is expected given that the prefit flux uncertainty includes beam pointing uncertainties, and this analysis is not expected to be able to strongly constrain them. Larger deviations of 5 mrad. were also tested, although are not shown, for which a strong bias is seen in the best fit χ 2 distribution.
Changes to the background predictions could also affect the result, although correlations between background and signal templates in the postfit covariances were weak (see the example in Figure 8), it may be expected that such changes will not have a large effect on the result. MINERvA have observed a deviation between Monte Carlo expectation and data as a function of reconstructed Q 2 (see Ref. [16] for a definition) for ν µ -CH andν µ -CH charged current quasielastic-like events. In particular, there is a strong suppression at low Q 2 values, which might correspond to the region of overlap between the signal and background templates. MINERvA has concluded, however, that the majority of this disagreement is due not to truly quasielastic events which would create background to a neutrino-electron scattering analysis, but rather to higher recoil processes such as pion production where the pion is observed in nuclear final state interactions [15,16]. In order to investigate this effect on the analysis, the input ν e -40 Ar andν e -40 Ar background events are modified from the nominal GENIE prediction by reweighting according to the observed MINERvA ratio, as a function of true Q 2 . The effect that the Q 2 modification has on the ν e background template is shown, as a ratio with the nominal template, in Figure 23, for both FHC and RHC modes. Figure 24 shows the effect on the weighted average flux values, with the MINERvA Q 2 distortion applied, for the nominal LAr detector in both FHC and RHC modes. In both cases, the bias on the weighted average flux distribution is very small except at high energies or in very low statistics bins.

F. Impact of radiative corrections
Although the signal ν-e − scattering process is well understood in principle, the radiative corrections applied in this analysis (described in Appendix A) are not included in GENIE, and may be refined by a more careful calculation later. To quantify the importance of the radiative corrections in the analysis, the study was repeated without radiative corrections applied, in the same manner as the previous studies into the effect of systematic uncertainties, although it is explicitly not a systematic uncertainty. Figure 25 shows the weighted average flux values produced with 10,000 throws of the input flux covariance, without radiative corrections applied (the default GE-NIE prediction), for the nominal LAr detector in FHC and RHC modes. The effect of removing the radiative corrections is to increase (decrease) the weighted average flux value in FHC (RHC) by 1-2%. Interestingly, the average χ 2 value from the 10,000 flux throws considered is slightly increased relative to the nominal case shown in Figure 16, the increase is larger than the systematic uncertainties considered in Section VI E, although not significantly. This can be understood because the absolute rate is well constrained by the analysis presented here, so a normalization-only change in the flux (e.g., a fully correlated shift between all energy bins) is more strongly constrained than some variations in the flux shape.

VII. INVERSE MUON DECAY
As previously noted in Section II, inverse muon decay (IMD), ν µ e − → ν e µ − has potential to constrain the high energy DUNE ν µ flux. The rate of such events in a 30 ton reference detector with a five year exposure is shown in Figure 26 as a function of both laboratory neutrino and muon energy. There are of order 6 × 10 3 events in such a sample produced with the neutrino mode beam, which would constrain the total rate of such high energy ν µ . The statistical sensitivity of such a sample to variations in this high energy ν µ spectrum is illustrated in Figure 27. The statistical uncertainty in such a sample would clearly allow either of two simulated neutrino spectral distortions, one of which increases the number of neutrinos above 15 GeV by 10%, and one of which increases the number of neutrinos above 24 GeV by 20%. However, distinguishing between the two scenarios would be more difficult.
In addition to the limited use of such a probe of the high energy ν µ flux for DUNE's primarily neutrino oscillation mission, the design of the near detectors will likely not be optimized to measure these spectra. Such high energy muons will not be contained in DUNE's near detectors, and measurement of the momentum by curvature in a magnetic field at such high muon energies will be difficult in a detector optimized for such measurements at significantly lower muon energies.

VIII. CONCLUSIONS
Because the neutrino-electron elastic scattering cross section is known, the flux can be extracted from the observed event rate. For realistic DUNE near detectors in the LBNF beam, it is possible to select a sample of many thousands of νe − → νe − events. In this work, we have investigated how well different potential DUNE near detector designs will be able to constrain the LBNF flux.
We found that given realistic mass constraints, a 30t liquid argon detector is able to perform better than 5t low density trackers, even ones with significantly better tracking resolution. This is due to its higher statistics, and despite the superior angular resolution of lower density detectors. This was also found to be the case even for a 5t detector with perfect electron reconstruction and background rejection.
With realistic systematic uncertainties, the uncertainty on the absolute neutrino flux in the 30t LAr detector is reduced from ∼8% to ∼2%. The uncertainty on the shape as a function of neutrino energy is also reduced by ∼20-30%. This is partially due to the fact that the flux shape is better known a priori than the absolute normalization. The improved reconstruction performance of high resolution detectors has a stronger impact on the flux shape constraint, as expected, but for realistic detector sizes, a large liquid argon detector still outperforms a smaller detector with better resolution. It seems that detector mass is the most important factor for making a ν-e − flux constraint, even for a 30t detector in the very intense LBNF beam. The intrinsic divergence of the beam is an important consideration which has the potential to limit the utility of a ν-e − flux constraint, and as such was included in this study.
As well as being able to reduce the neutrino flux uncertainties, we demonstrated that a ν-e − sample with a large liquid argon detector is capable of identifying a large variation in the neutrino spectrum outside of predicted uncertainties, and is not biased to the input assumptions about the flux, despite the inability to directly distinguish neutrino flavor with a ν-e − sample.
fee Shop, Qui Nhon, Vietnam, for their gracious hospitality, their dehumification devices that kept the air no worse than unpleasantly tepid, and their limitless supply of tasty iced beverages. We also acknowledge the valuable feedback from many members of the DUNE collaboration as this work was presented in the context of the near detector design studies.
In the second case, where y ≡ (T e + E γ )/E ν , the mod-ifications to the cross-section are more straightforward: − 2y 9(1 − y) 2 + 23 72(1 − y) 2 (A8) X 3 is not available in the calculation of Ref. [30], although it has been recently calculated by the authors of Ref. [11]. However, since X 3 enters into Eqns A1 and A2 only in terms multiplied by m e /E ν , it can be safely neglected.
The direction of the electron, however, in any detector under consideration is most likely to be measured as the electron's direction. All of these calculations are done assuming collinear emission of photons along the lepton angle. Within that approximation, where ≡ E γ /(T e + E γ ). Note that corrections from photon emission occurs only multiplied by the small θ 2 e and is thus negligible. This implies that there is no significant effect due to real photon radiation on the reconstructed neutrino energy inferred from the electron angle and clustered electron plus photon energy, such as in a liquid argon TPC or other calorimetric detector.