Measurement of the Double-Differential Muon-neutrino Charged-Current Inclusive Cross Section in the NOvA Near Detector

We report cross-section measurements of the final-state muon kinematics for \numu charged-current interactions in the NOvA near detector using an accumulated 8.09$\times10^{20}$ protons-on-target (POT) in the NuMI beam. We present the results as a double-differential cross section in the observed outgoing muon energy and angle, as well as single-differential cross sections in the derived neutrino energy, $E_\nu$, and square of the four-momentum transfer, $Q^2$. We compare the results to inclusive cross-section predictions from various neutrino event generators via $\chi^2$ calculations using a covariance matrix that accounts for bin-to-bin correlations of systematic uncertainties. These comparisons show a clear discrepancy between the data and each of the tested predictions at forward muon angle and low $Q^2$, indicating a missing suppression of the cross section in current neutrino-nucleus scattering models.


INTRODUCTION
Neutrino scattering on nuclei is a rich topic with many challenges, both experimentally and theoretically.Experimentally it is challenging to produce a wellcharacterized source of neutrinos and to collect high statistics with high-resolution detectors.Accordingly, many recent inclusive neutrino-nucleus scattering measurements are limited by large statistical and/or systematic uncertainties [1][2][3][4].Theoretical challenges arise from a lack of accurate models that are valid across a large range of energies and account for the initial state of the nuclear environment and final-state interactions [5].
In current and future long-baseline neutrino flavor oscillation experiments, beams of muon (anti)neutrinos are used to precisely measure the rate of muon (anti)neutrino disappearance and electron (anti)neutrino appearance [6][7][8][9].Weak charged-current (CC) interactions, in which a charged lepton is produced in the final state, are used to identify neutrino flavor and measure the neutrino energy.The accuracy of these measurements depends explicitly on the kinematics of the lepton and hadrons visible in the final state.
To relate these final state observables to the energy of the neutrino, accurate knowledge of neutrino-nucleon interaction cross sections and the dynamics of the propagation of particles through nuclear matter is necessary.Many of the uncertainties in neutrino oscillation parameters that arise from limited understanding of neutrino cross sections are reduced by using a two detector scheme with a near detector placed close to the beam source to characterize interactions prior to oscillation [10,11], and a far detector placed much farther away to measure the oscillated neutrino spectra.However, the near and far detectors are typically substantially different in size and have differing acceptances of the final-state particles produced by neutrino interactions in the detector.Knowledge of kinematic distributions of the final-state leptons is crucial to correctly account for differences in event selection efficiency and purity between the two detectors.In practice, experiments rely on neutrino event generators for this knowledge.
Neutrino interactions are typically characterized by the type of target (e.g., individual nucleons, pairs of nucleons, the nucleus as a whole or electrons) and the particles produced in the interaction.At around 1 GeV in neutrino energy, quasielastic (QE) scattering dominates, in which the neutrino scatters off a single nucleon, producing a lepton and a single unbound nucleon.At these energies, meson exchange currents (MEC) between pairs of correlated nucleons resulting in 2-particle-2-hole (2p2h) interactions also significantly contribute to the neutrino scattering rate [12][13][14].Around 2 GeV, resonant (Res) interactions contribute significantly to the total cross section.In these interactions, intermediate hadronic excited states are created inside the nucleus (predominantly those associated with ∆(1232) resonances) that decay to a baryon and a meson.At energies above 3 GeV, shallow-and deepinelastic scattering become more prevalent.Other rarer interactions, such as neutrino-electron and coherent scattering off the entire nucleus (COH) also contribute to the total cross section.
Nuclear effects play a significant role in the initial and final states of the interaction.In addition to introducing new primary processes such as MEC as noted above, at GeV energies the initial state of the nucleus influences the kinematics of the particles produced in the interaction [13][14][15].Furthermore, these particles must traverse the nuclear medium, during which scattering and interactions may occur.These final-state interactions (FSI) alter the kinematics and possibly the composition of the final state [5].
The inclusive cross section, σ incl , is the sum of the cross sections of all of the individual processes.As such, predictions of the inclusive cross section must properly combine the individual processes, including interference terms.Measurements of inclusive cross sections serve to constrain the quantum-mechanical sum of these processes, as well as their dependence on the neutrino energy (E ν ) and square of the four-momentum transfer from the lepton system (Q 2 ), and the impact of final-state interactions.
We report the flux-integrated double-differential inclusive cross section of neutrino-nucleus CC interactions in the NOvA near detector, ν µ + A → µ + X, where A is a target nucleus (see Table I) and X represents all other final state particles.The measurement is differential with respect to the final-state muon's kinetic energy and angle relative to the neutrino beam direction.We also report the inclusive cross section as a function of the derived E ν and Q 2 , integrated over the range of muon kinematics reported in the double-differential measurement.

THE NOVA EXPERIMENT
NOvA is a long-baseline neutrino experiment [6] designed to measure neutrino flavor oscillations.A 96% pure muon-neutrino beam is produced at Fermilab.Two functionally identical detectors are directly exposed to the beam: the near detector located 1 km downstream of the beam target, and the far detector located 810 km away from the target near Ash River, Minnesota.The primary measurements of electron (anti)neutrino appearance and muon (anti)neutrino disappearance provide constraints on the neutrino mixing parameters, θ 23 and ∆m 2  32 , the neutrino mass ordering, and the CP-violating parameter, δ.The high statistics neutrino and antineu- trino samples gathered at the near detector constrain the flux and neutrino cross-section parameters for the oscillation analyses, and are also ideal for measurements of various neutrino interaction cross sections.
Neutrinos for NOvA are provided by the Fermilab NuMI beam [16].The Fermilab Main Injector protons at 120 GeV strike a graphite target, producing pions and kaons.These hadrons are focused by two magnetic horns and directed towards a 650 m drift region where they decay to produce primarily muons and muon neutrinos.The horn polarity can be changed to focus positive (negative) mesons and produce a primarily (anti)neutrino beam.The NOvA detectors are located 14.6 mrad offaxis from the central beam direction, resulting in an incident neutrino energy spectrum narrowly peaked at 1.8 GeV. Figure 1 shows the flux at the NOvA near detector in the neutrino beam configuration.The neutrino beam includes a 1.8% intrinsic νµ component coming from opposite-sign meson decay in the energy range of interest for this measurement, between 1 and 5 GeV.There is also an electron neutrino and antineutrino contribution of 0.7% in this energy range.
The NOvA near detector is a tracking calorimeter with 193 t of active mass, located 100 m underground.The detector is composed of planes of hollow cells made from a custom formulation of extruded PVC.The planes are segmented into 3.9 cm wide cells that are 3.9 m long; the depth of each plane in the beam direction is 6.6 cm.The planes are alternated in horizontal and vertical orientations perpendicular to the beam, allowing full 3D tracking for 12.7 m along the beam axis.Each cell is filled with liquid scintillator, a blend of 95% mineral oil and 5% pseudocumene with trace concentrations of wavelength shifting fluors.The resulting composition by mass is about 63% scintillator and 37% PVC with nuclear targets for neutrino interactions in the detector as described The downstream end of the detector is a "muon catcher" designed to improve containment of muons produced in neutrino interactions up to ∼2.5 GeV.The muon catcher consists of 10 layers of 10 cm thick steel absorbers interleaved with pairs of PVC/scintillator.The muon catcher planes span the full detector width and the lower 2/3 rds of the detector height.

SIMULATION
Simulation is used in this analysis to calculate the integrated flux, selection efficiencies and purities, estimate energies, and effects of detector resolutions.The analysis also relies on simulation to optimize event selection criteria and assess various systematic uncertainties that can impact event rates and selection efficiency and purity.The simulation is a chain of steps that begins with the generation of the neutrino beam and transport of all particles through the beamline to the detector.Interactions of the neutrinos with the detector are then generated, after which the final-state particles are transported through the detector.The generation, detection and digitization of light in the detector are the final steps of the simulation chain.Each step of the simulation chain, described below, is matched to the real data-taking conditions in beam intensity and total protons-on-target, wherever appropriate.
The NuMI flux predictions start with a detailed simulation of the beamline components and the hadronic showers that follow the primary proton striking a long graphite target until the mesons decay to neutrinos.The simulation is based on GEANT4 v9.2.p03 [17] with the FTFP BERT hadronic model.The hadron production model is adjusted using the PPFX package, which uses external measurements on thin targets with the procedure outlined in Ref. [18].The NuMI flux prediction for The ratio in the lower panel is the size of the flux correction using PPFX [18] with respect to the GEANT4 v9.2.p03 [17] FTFP BERT hadronic model.
the neutrino beam mode at the NOvA near detector is shown in Fig. 2.
The simulated neutrino flux is passed through a detailed description of the NOvA near detector geometry, including surrounding rock, where interactions are simulated with the GENIE v2.12.2 [19,20] event generator.The initial state is simulated via the default Smith and Moniz global relativistic Fermi gas (RFG) model [21].Short-range nuclear correlations in the initial state [22] are accounted for by the addition of a high-momentum tail of the Fermi momentum distribution for single nucleons [23].The QE interactions are simulated according to the formalism of Llewellyn Smith [24].2p2h interactions are simulated using the Empirical MEC model [25].Charged-current Res interactions are simulated via the Rein and Sehgal model [26].Inelastic scattering over a large range of hadronic invariant masses, resulting in a range of final state hadrons, is simulated using the Bodek-Yang scaling formalism [27] coupled to a custom hadronization model [28] and PYTHIA 6 [29].Chargedcurret COH interactions are simulated using [30,31].
The GENIE output has been adjusted to incorporate advances in theory and experimental data [12].These modifications include adjustments to the CCQE and nonresonant pion production interactions based on reevaluated bubble chamber measurements; improved nuclear models of CCQE kinematics; and suppression at low Q 2 of resonant pion production.After applying these modifications, differences in the shapes of the energy and threemomentum transfer (q 0 ,| q|) distributions of near detector data and simulation were used to tune the Empirical MEC model.These adjustments significantly enhance the agreement between selected muon-neutrino candidates in the NOvA near detector data and simulation across multiple kinematic variables such as visi-ble hadronic energy and reconstructed four-momentum transfer squared.GEANT4 v10.1.p3is used to simulate energy deposited in the NOvA near detector from the particles generated by neutrino interactions.A custom simulation tuned to reproduce measured scintillator response and fiber attenuation properties is then used to model and transfer scintillation and Cherenkov light [32].Test-stand measurements are used to tune the Birk's suppression of the scintillation light and to validate the response of the readout electronics in the simulation [33].

EVENT RECONSTRUCTION AND CALIBRATION
Energy deposits (hits) in the detector are recorded with pulse height, time and channel location information.Cell-to-cell variations in pulse height are first corrected using through-going muons, followed by the calibration of absolute energy deposition using minimumionizing portions of stopping cosmic ray muon tracks.The reconstruction of neutrino interactions first clusters hits that are correlated in space and time [34].These clusters of hits are all assumed to be associated with a single neutrino interaction, referred to as an event.Hits in an event are then grouped into possible particle trajectories (tracks) via a Kalman filter-based algorithm in both the horizontal and vertical two-dimensional detector views [35].Three-dimensional tracks are formed by combining tracks from the two views based on their overlap in the longitudinal direction.The track reconstruction algorithm assumes a start point at the most upstream hit, and requires a minimum of 4 hits in each detector view.A different algorithm [36] is used to form particle trajectories (prongs) from hits associated with a reconstructed vertex, and requires a minimum of 1 hit in each detector view.As described later, tracks and prongs are used for different purposes in this analysis.

EVENT SELECTION
Candidate events are required to have a reconstructed track which crosses more than 4 contiguous planes, made up of hits in at least 20 unique cells.The candidate muon track, described below, is required to start inside a 2.7 m × 2.7 m × 9 m fiducial volume inside the detector upstream of the muon catcher.We require all tracks and prongs identified in the event to have stopped several cm before reaching any detector edge to ensure containment of all the neutrino energy.We further require that no track or shower other than the selected muon enter the muon catcher.These criteria select 23% of signal events.
The signature of a ν µ CC interaction is the presence of a muon in the final state.This analysis implements a multivariate muon identification algorithm, MuonID, based on energy deposition and scattering observables (see Fig. 3).For energy deposition, we use the difference between log-likelihood functions based on the dE/dx of a muon and a pion, the average dE/dx in the cells of the last 10 cm of the track trajectory, and average dE/dx in the cells of the last 40 cm of the track trajectory.We also use distributions of the difference between log-likelihood functions based on the angular deflections along the trajectory of the reconstructed track for muons and pions.These reconstructed variables are used as input to a boosted decision tree (BDT) algorithm, the output of which is a MuonID score.The BDT is trained on all true muon tracks and true non-muon tracks using reconstructed simulated neutrino interactions that have passed the preselection criteria described above.The samples used to train and test these algorithms are drawn from non-overlaping subsamples that each comprise 10% of the overall simulated sample.The distributions of the highest MuonID score in signal and background events passing the preselection are shown in the left-side plot of Fig. 4.
As this measurement is systematically limited, we optimize the MuonID selection criteria by minimizing a figure-of-merit (FOM) that is approximately the fractional uncertainty on the total cross section: where is the selection efficiency and P is the selection purity.The sources of uncertainty considered for the selection criteria optimization are neutrino interaction modeling, energy scale uncertainties, and the modeling of light generation and propagation in the detector.These sources of systematic uncertainty, described in more detail in Sec. 8, have the greatest impact on muon identification.The right plot in Fig. 4 shows the uncertainty on the purity and efficiency, and the FOM as a function of the minimum MuonID value in the event.Signal events with MuonID greater than 0.24 are retained as candidate ν µ CC interactions, resulting in an overall 98% selection efficiency and overall 97% selection purity after the previously described selection criteria are applied.The muon is correctly identified in 98.8% of signal events.The neutrino interaction vertex is taken as the most upstream position of the selected muon track.

ENERGY RECONSTRUCTION, RESOLUTION AND BINNING
The muon and muon-neutrino energy estimators developed for this analysis rely on the simulation to relate the muon energy to the length of the reconstructed muon track.Muons that stop before reaching the muon catcher are reconstructed with a typical energy resolution of 4%; those that stop in the muon catcher have a resolution of 5-6%.
We reconstruct the visible hadronic energy as the sum of calibrated energy of hits in the event that are not as-  sociated with the muon track, plus any additional energy that may be deposited by hadrons on and near the start of the muon track.The latter is reconstructed by subtracting the energy of a minimum ionizing particle from the first few planes of the muon track.We then use the simulation to convert the visible hadronic energy to an estimate of E avail [13], the total true energy of the visible hadrons in the final state.
The cross section is reported as a function of the directly observed kinetic energy of the muon, T µ , and the cosine of the angle of the muon with respect to the neutrino beam direction, cos θ µ .The cross section is also reported as a function of model-dependent E ν and Q 2 .The combination of the muon kinematics and E avail is mapped to E ν using simulation, and then the combination of the reconstructed E ν and muon kinematics are mapped to Q 2 using simulation.All bins are at least as wide as the resolution estimated from simulated signal events that pass the selection.The average T µ resolution is 50 MeV, and the resolution of the average muon angle is typically less than 4 • .We use 20 equal-sized bins from 0.5 GeV to 2.5 GeV for reconstructed T µ , and 13 variable-sized bins for reconstructed cos θ µ between 0.5 and 1.The choice of the variable-sized binning in cos θ µ accounts for both resolution and statistics, with smaller bins in the most forward angles.The binning can be seen in Fig. 5, discussed below.

THE MEASUREMENT AND RESULTS
The results presented in this paper use 8.09 × 10 20 protons-on-target (POT) collected between August 2014 and February 2017 in the neutrino beam configuration.The double-differential cross section is determined as: Here N sel is the number of selected events, P is the selection purity (the fraction of signal events among selected events), is the selection efficiency (the fraction of signal events selected), φ is the integrated neutrino flux, N target is the number of nucleon targets in the fiducial volume, ∆ cos θ µ is the width of the angle bin, and ∆T µ is the width of the muon kinetic energy bin.An unfolding matrix, U −1 ij , is used to relate the reconstructed observable in bin j to the true observable in bin i.As seen in Fig. 5, the analysis receives non-negligible contributions from different interaction modes, each with varying amounts of hadronic energy in the final state.Hadrons in the final state can influence the purity, unfolding and efficiency, for example when charged pions are misidentified as muons, or if a hadronic shower hides the presence of the muon, or if the hadronic system is too close to the edge of the detector and the event fails the containment criteria.Therefore, a three-dimensional space involving the muon kinematics and E avail is used in applying the purity, unfolding, and efficiency corrections to reduce potential model dependences on the finalstate hadronic system.Ten 250 MeV-wide bins (and one overflow bin) are used for E avail .The corrected threedimensional result is then integrated over E avail .
In order to compare this measurement to theoretically based predictions, we use the D'Agostini iterative unfolding algorithm [37] as implemented by the RooUnfold package [38] to correct for bin-to-bin migrations of selected signal events.The number of iterations performed is a regularization parameter that serves to reduce extreme variations in the unfolded distribution that are consistent with the data given the predicted smearing but are implausible given the underlying physics of the system.We choose the number of iterations which minimizes the weighted mean of the relative bias and variance across all bins using independent simulation samples with randomized systematic shifts.The generation of systematically shifted simulations is described in Sec. 8.The nominal simulation prediction for the response matrix was used to unfold the shifted simulation data sets.We found the optimum number of iterations for this analysis to be between 2 and 4 for a variety of shifted simulation data sets where both signal and background normalizations and shapes were systematically shifted, and chose 3 iterations for the unfolding applied to the data.
The purity, P , and the efficiency, , are shown in Fig. 6 vs. T µ for each bin of cos θ µ .Curves are drawn separately for various representative ranges of E avail .Each angular grouping shows a clear dependence of the purity on the muon kinetic energy.At low T µ our selection suffers from contamination by NC interactions.This effect is more evident at higher available energy, as a higher fraction of hadronic activity increases the chances of misidentification.The efficiency increases with increasing cos θ µ , as at larger angles the muon is more likely to escape via the side of the detector or less likely to be clearly separated from hadronic activity in the detector and therefore less likely to be reconstructed as a track or identified as a muon.The efficiency decreases as a function of muon kinetic energy as higher energy muons are less likely to be contained in the detector.We also observe a clear dependence on E avail , as a larger fraction of hadronic activity makes the event reconstruction and identification of the muon more difficult.Comparisons of the purity and effi- ciency of the event selection with and without the NOvA tune of the simulation were found to be in agreement within systematic uncertainties.
Figure 7 shows the extracted double-differential cross section in slices of muon angle.Figure 8 shows the extracted single differential cross section vs. Q 2 and vs. of E ν , restricted to the phase-space of the doubledifferential measurement.Efficiency, purity, and unfolding corrections are simple functions of Q 2 (E ν ) for these derived quantities.The data are presented with total and statistical error bars in the plot, and the values are also available in the table in Appendix A, and in electronic format at the NOvA Experiment Data Releases webpage [39].The data are compared to predictions from GENIE v2.12.2 with and without the tune described in Sec. 3. We observe better than 5% agreement between the data and the NOvA-tuned GENIE v2.12.2 prediction across all muon angle slices, although small discrepancies are still present, especially at low (∼1 GeV) muon kinetic energies and very forward angles.The tuning procedure does not significantly impact the predictions at larger muon angles and so the untuned predictions are very similar to the tuned predictions.However a clear discrep-ancy between the data and the untuned GENIE v2.12.2 prediction is evident in the three most forward-going angle bins (cos θ µ > 0.96).As shown in Fig. 5, these three bins are heavily dominated by QE, MEC and resonant interactions and are particularly sensitive to the low-Q 2 suppression discussed in Sec. 3. The Q 2 discrepancy between the data and the GENIE v2.12.2 predictions is shown in the left side of Fig. 8, where the data imply a need for an even stronger suppression of the cross section at very low values of Q 2 than is currently achieved via the simulation tuning procedure.The Q 2 discrepancy is washed out as a function of neutrino energy, so there is broad agreement between the data and predictions in the right-side of Fig. 8.

UNCERTAINTIES
Several sources of systematic uncertainty impact this measurement: the neutrino flux prediction, detector response, muon energy scale, muon angle, normalization, modeling of neutrino-nucleus interactions, and modeling of neutron interactions in the detector.In general, for each source of uncertainty, we use the difference between our nominal simulation and systematically modified simulations to estimate the uncertainty on the selection efficiency and purity.If the source of systematic uncertainty could impact the reconstruction or particle identification algorithms, then the effect is applied to the same simulated neutrino interactions at the relevant point in the simulation-reconstruction chain, and the effect is propagated through the reconstruction and analysis chain.
Otherwise, the impact of a systematic source is estimated by applying weights to events in the simulation.In all cases, the migration matrices used in the unfolding procedure are recalculated for each systematic variation.
For each individual source of systematic uncertainty, a systematically shifted "universe" is simulated with a ±1σ shift to the systematic source parameter.Calibration and muon energy scale are examples of systematic uncertainties for which this approach is used.Other uncertainties, such as neutrino cross-section modeling and flux, are impacted by many sources and are calculated with a multi-universe method.In this method, a hundred or more universes are generated where parameters influencing the uncertainty are drawn from a normal distribution, with a width that corresponds to the 1σ uncertainty on each systematic source parameter.
Uncertainties in the neutrino flux prediction arise from the modeling of hadron production in the target, horns and decay pipe, and from the modeling of the beam optics.The hadron production uncertainty on the neutrino flux after the adjustments explained in Sec. 3 is ∼7% at the spectrum peak.This uncertainty is dominated by interactions for which there are no relevant external data to be included in the adjustment procedure (mostly meson and proton elastic and quasielastic scattering).Uncertainties in beam optics are incorporated by propagating uncertainties in the alignment and focusing of beamline elements; this uncertainty is ∼4% at the peak.
Detector response uncertainties include uncertainties in the calibration of the visible hadronic energy scale and simulation of light production and transport from the liquid scintillator and wavelength-shifting fibers to the photodetectors.A 5% difference in the recorded energy deposition as a function of distance traveled of candidate proton prongs measured between simulation and data is used as the uncertainty in the hadronic energy response.We use systematically shifted simulation samples where the absolute energy scale is shifted by ±5% to evaluate the impact on this analysis.An observed non-uniformity in the calibrated energy response as a function of distance from the readout is included as a calibration shape uncertainty.The uncertainty on the light model arises from the uncertainty on overall light yield of the scintillator and the efficiency with which Cherenkov photons are absorbed by the scintillator and re-emitted at wavelengths that can be detected.A simulation sample where Cherenkov light production is disabled is used to assess an upper limit on the uncertainty on this aspect of the light model.Uncertainties in the muon energy scale arise from modeling the energy loss of muons in the detector.A detailed analysis of muon energy loss in the NOvA near detector material composition in Geant4 indicates a ±0.8% uncertainty for the portion of the track that traverses the fully active region, and ±1.2% for the portion of the track that traverses the muon catcher [40].We conservatively assume the worst-case scenario and scale the reconstructed muon energy in these fractions by either all positive or all negative directions in both regions of the detector in assessing this uncertainty.
Uncertainties in the muon angle arise from misalignments of the PVC cells in the near detector.To estimate the impact of these misalignments on the muon direction, an alternative simulation sample was generated with randomly shifted cell positions according to the construction tolerances of the detector.A comparison of this systematically shifted simulation to the nominal simulation shows a 2.5 mrad spread in the reconstructed muon angle distribution, and negligible spreads in the muon and hadron energy distributions.We implement a 2.5 mrad systematic shift to the reconstructed angle of the muon to determine the impact on the measured cross section.
Normalization uncertainties in the measured cross section arise from uncertainties in the detector mass, integrated POT exposure and modeling of beam intensity.Data from the manufacturing and construction processes are used to constrain the uncertainty on the mass of the detector to 0.28%, and the uncertainty on the POT accrual in the NuMI beamline is 0.5% based on measurements of beam current through a toroid magnet.The simulation accounts for time-dependent variations of beam intensity.An observed ∼2% difference between shapes and normalizations of data and simulation selection efficiencies as a function interaction vertex position in this analysis is used as the uncertainty due to beam intensity modeling effects on the normalization.The combined uncertainty on the normalization of the reported cross section is 2.1%.
We use a reweighting approach to estimate the im- Figure 7. Extracted double-differential cross section, in slices of muon angle.The data are presented showing statistical and total uncertainties, and compared to GENIE v2.12.2 -NOvA Tune [12] (solid red line) and GENIE v2.12.2 -Untuned [19,20] (dashed red line).The inner error bars are from statistics, and may not be visible on this scale.).The data are presented showing statistical and total uncertainties, and compared to GENIE v2.12.2 -NOvA Tune [12] (solid red line) and GENIE v2.12.2 -Untuned [19,20] (dashed red line).
The inner error bars are from statistics, and may not be visible on this scale.
pact of neutrino-nucleus scattering uncertainties.The weights applied are a mix of NOvA-specific uncertainties and uncertainties available from the GENIE event generator [20].The NOvA-specific uncertainties include a 5% uncertainty on the value of the CCQE M A parameter and a 100% one-sided uncertainty on the Q 2 suppression of resonant pion production applied to the simulation.For MEC interactions, uncertainties in the fraction of target nucleon pairs (np vs. pp) in the nucleus and the dependence of the MEC cross section as a function of q 0 and q 3 are taken into account.Additional NOvA-specific uncertainties are included for DIS interactions.Further details of the NOvA-specific uncertainties are described in [12].An energy uncertainty is assigned to the detector's response to neutrons.This uncertainty is driven by comparison of data to simulation in a neutron-rich subsample of the antineutrino dataset.An excess of neutrons with low visible energy is observed.A sample where one third of the neutron candidates with energy below 40 MeV had their visible energy scaled down by a factor of 3.6 produces better data-simulation agreement.The sample is used to set the size of a conservative two-sided neutron response uncertainty.
Bin-to-bin correlations from all sources of systematic uncertainty are derived from using the difference between 10 5 systematically-shifted simulations and the nominal simulation to calculate a total systematic uncertainty covariance matrix.The unfolding procedure also induces small bin-to-bin correlations.We calculate the statistical covariance matrix using a multi-universe procedure similar to that described above, with 4000 toy simulations with Poisson-fluctuated event counts in each measurement bin in reconstructed space.The total uncertainty covariance matrix is a linear sum of the total systematic and the statistical covariance matrices.Table II shows the breakdown of the weighted average fractional uncertainties and correlations in the double-differential crosssection measurement.The weighted average fractional uncertainty is defined as where i is a measurement bin, V is the covariance matrix, and σ is the measured double-differential cross section.The weighted average correlation is defined as where i is a measurement bin, j is a different measurement bin (so that diagonal elements are excluded), and C is the correlation matrix.The dominant source of uncertainty comes from the flux prediction.As the average correlation indicates, this is mainly a normalization uncertainty and is significantly reduced in the shape-only analysis, where normalization differences have been removed.We note that the statistical uncertainties are at the level of a few percent per bin, and that the interaction modeling uncertainties are subdominant.The normalization uncertainties that are 100% correlated across all bins are removed in a shapeonly covariance matrix.The typical total uncertainty of our measurement is around 12%, which is reduced to 7% in the shape-only analysis.

COMPARISONS TO GENERATORS
The results of the double-differential and single differential cross-section measurement are presented in this section and compared to GENIE versions v2.12.2 and v3.00.06,NEUT v5.4.0 [41], NuWro 2019 [42,43] and GiBUU 2019 [44,45].Table III lists the models used in each for the initial state, interaction modes, and final state interactions in the generators.GENIE v2.12.2 is the neutrino event generator used in the simulation for this analysis and is described above in Sec. 3. GENIE v3.00.06 is a more recent version of GENIE, and we use a configuration chosen by the NOvA experiment for its 2020 oscillation analysis, N18_10j_02_11a, a combination of G18_10j_00_000 and G18_10b_02_11a which in practice results in predictions that are nearly identical to the out-of-the-box predictions from the G18_10b_02_11a tune.We note however that the GENIE v3.00.06 tune shown here has no other NOvA-specific tuning applied.Whereas all other event generators use a local Fermi-gas model (LFG) for the initial state, GENIE v2.12.2 uses a relativistic Fermi-gas (RFG) model.QE and MEC interactions are implemented via the València group model from Nieves, et al. [46,47] in GENIE v3.00.06 and NEUT v5.4.0.NuWro implements QE interactions based on the Llewellyn Smith [24] model with an additional Random Phase Approximation (RPA) suppression, but implements MEC interactions based on the València model.Resonant interactions are based on the Berger-Sehgal [48] model, and DIS interactions use PYTHIA 6 [29] in GE-NIE v3.00.06 and NuWro 2019 and PYTHIA 5 in NEUT v5.4.0.GiBUU implements its own unique model for neutrino interactions across the QE and Resonant regions based on many of the same principles as the models mentioned above [49,50].Final-state interactions are implemented via a variety of models, including that of Oset et al. [51], cascade models in GENIE and the Boltzmann-Uehling-Uhlenbeck (BUU) equations in GiBUU.It is worth noting that most of the models listed above use form factors and cross sections extracted from very similar data sets, and in principle should be highly correlated.However, as we show below, and has been noted elsewhere (see e.g., [5]), the inclusive charged-current neutrino-nucleus cross section predictions from these different generators differ considerably, likely due to differences in implementation.
The top panel of Fig. 9 shows the ratio of the extracted double-differential cross section to the GENIE v2.12.2 -NOvA Tune [12] prediction, in slices of muon angle.The outer error bars of the data represent total   uncertainties, while the inner error bars of the data are statistical only.The solid histograms are ratios of the predictions from different neutrino event generators to the GENIE v2.12.2 -NOvA Tune prediction.In the lower panel of Fig. 9 the predictions are first area-normalized to the data across the reported double-differential measurement space before taking the ratio with respect to the GENIE v2.12.2 -NOvA Tune prediction, and the outer error bars of the data represent shape-only uncertainties.These comparisons indicate 5-10% agreement between the measurement and the various generators at high-angle slices.Discrepancies become more apparent at more forward-going angles and lower muon energies.
Figure 10 shows similar comparisons of the differential cross section as a function of Q 2 and the cross section as a function of E ν .These model-dependent variables are calculated only in the muon kinematics space specified by the double-differential measurement.The top plots show unmodified predictions and the data with total error bars.The bottom plots show predictions that are area-normalized to the data and the data with shape-only error bars, where the normalization uncertainties that are 100% correlated across bins have been subtracted.As was the case for GENIE v2.12.2, the large suppression at Q 2 < 0.1 GeV 2 is not described by any of the generators.As the generators use very similar models for interactions that contribute most at low values of Q 2 (QE and MEC), this strongly indicates that some additional suppression of the cross section at low Q 2 is lacking from the underlying theory.Furthermore, many of the predictions prefer a stronger suppression of the cross section at Q 2 values ranging from 100-800 MeV than is observed in the data, which suggests the need for improved modeling of resonant interactions.The differences between the data and predictions as a function of Q 2 are washed out when looking at the cross section as a function of the neutrino energy, and we see good agreement between our measurement and most neutrino generators.
In order to make a more quantitative assessment of the agreement between our measurement and the various event generators, we calculate the global χ 2 between our measurement and different generators across all measurement bins.We use the systematic uncertainty covariance matrix described in Sec. 8 to account for binto-bin correlations in the χ 2 calculation.Table IV shows a summary of the global χ 2 calculations for both total and shape-only uncertainties.The normalization factor used to area-normalize the predictions to the data for the shape-only comparisons is also shown.GENIE v2.12.2 with the NOvA tune results in the best χ 2 , however we note that the χ 2 per degree of freedom (dof) of ∼2 (there are 158 degrees of freedom) is yet another reflection of the remaining discrepancies between the measured and tuned predicted cross section.As expected from the data-generator comparisons in Fig. 9, the global χ 2 s are very high, but vary across generator predictions.The global χ 2 s for the shape-only comparison are even larger, which implies that shape differences are significant and that a simple normalization correction to the predictions is insufficient to reduce the discrepancies.The combination of the differences seen in Fig. 9 with the information on the interaction types in Fig. 5 emphasizes the regions of muon kinematic phase space and perhaps the particular models in each generator that need the most attention by the neutrino-nucleus scattering community.

CONCLUSION
We have presented a measurement of the doubledifferential ν µ CC inclusive cross section in the NOvA near detector in 158 bins of muon momentum and angle.The measurement applies purity, unfolding and efficiency corrections based on muon energy, muon angle and amount of observable hadronic energy in the detector, reducing the neutrino-nucleus interaction model dependence on the measurement.The measured cross sections and the covariance matrices are available in digital format on the NOvA Experiment Data Releases webpage [52].The weighted average fractional total uncertainty of 12% is driven primarily by a 9.1% flux normalization uncertainty.The flux normalization uncertainty is expected to decrease by more than a factor of two over the next few years as new constraints become available from external hadron production experiments such as NA61/SHINE (eg, [53,54]) and EMPHATIC (eg, [55,56]) as well as neutrino-electron scattering measurements in the NOvA near detector.The weighted average fractional shape-only uncertainty of 8.1% is driven by muon and hadronic energy scale uncertainties.Comparisons to generator predictions are made by calculating χ 2 using both the total and shape-only covariance matrices.There is an apparent tension between the measurement and predictions at very forward angles, consistent with a large observed discrepancy between the measurement and predictions at very low Q 2 .This discrepancy at low Q 2 is seen in all generators regardless of normalization uncertainties, and is consistent across all neutrino event generators to which the data are compared.The region of phase space covered by very forward muon-scattering angles receives contributions from QE-like and resonant interactions.Consequently, the data strongly suggest that a fundamental component responsible for greater low-Q 2 suppression of the cross section is missing from the interaction models.We note too that since a presentation of preliminary results of this analysis [57], the Giessen group has modified the resonant and shallow-inelastic scattering region in the GiBUU simulation improving agreement with our data [58], exemplifying the rapid pace of neutrino generator development and the need for additional data.Future measurements of neutrino interactions by the NOvA collaboration (e.g.see [59,60]) aim to isolate the exclusive final states that could be contributing to the large discrepancies observed in the inclusive channel presented in this paper.
Table IV.Summary of global χ 2 calculations for different neutrino generators for the double-differential cross-section measurement across 158 bins of muon energy and angle.The scale factors and χ 2 s for shape-only comparisons are also shown.The χ 2 calculation accounts for bin-to-bin correlations using the statistical and systematic covariance matrix described in Sec. 8

Figure 2 .
Figure 2. Muon-neutrino flux spectrum at the NOvA near detector below 5 GeV.The shaded band represents the total 1-σ flux uncertainty (hadron production and beam optics).The ratio in the lower panel is the size of the flux correction using PPFX[18] with respect to the GEANT4 v9.2.p03[17] FTFP BERT hadronic model.

Figure 3 .
Figure 3. Simulated muon (hashed blue) and non-muon (dashed red) track distributions in: dE/dx log-likelihood differences between that of a muon and a pion (top left), multiple scattering log-likelihood differences (top right), average dE/dx in last 10 cm (bottom left) and average dE/dx in last 40 cm (bottom right) used in the MuonID selector.Muon distributions are normalized to data exposure (8.09 × 10 20 POT), non-muon distributions are normalized by area to the muon distributions.

Figure 4 .
Figure 4. Left: Stacked distributions of the simulated maximum MuonID in each reconstructed event for signal (hashed blue), neutral current (NC, dashed red) and electron-neutrino (solid green) backgrounds.Right: Fractional uncertainty on the selection efficiency, selection purity and FOM versus the required largest MuonID value per event.Candidate νµ CC events are retained with a requirement of largest MuonID value in an event > 0.24.

Figure 5 .
Figure 5. Relative cross-section contributions for different interaction modes (QE -purple dotted filling, MEC -red right diagonal lines, Res -orange left diagonal lines, DIS -green horizontal lines, Other -all non νµ CC contributions in blue vertical lines) in the NOvA-tuned version of GENIE v2.12.2 as a function of Tµ for each bin of cos θµ.

Figure 6 .
Figure 6.Event selection purity (top) and efficiency (bottom) versus Tµ for each bin of cos θµ in representative ranges of E avail .Solid lines are for 0.00 GeV < E avail < 0.25 GeV, dashed oranges line are for 1.00 GeV < E avail < 1.25 GeV, dotted green lines are for 1.50 GeV < E avail < 1.75 GeV, and dashed-dotted blue lines are for 2.50 GeV < E avail < 120 GeV.

Figure 8 .
Figure 8. Single-differential cross section as a function of Q 2 , (left) and Eν (right) calculated only in the muon kinematics space specified by the double-differential measurement (see Sec. 6).The data are presented showing statistical and total uncertainties, and compared to GENIE v2.12.2 -NOvA Tune[12] (solid red line) and GENIE v2.12.2 -Untuned[19,20] (dashed red line).The inner error bars are from statistics, and may not be visible on this scale.

Figure 9 .
Figure 9. Extracted double-differential cross section divided by the GENIE v2.12.2 -NOvA Tune prediction.Ratios are shown in slices of muon angle and compared to the ratio obtained from GENIE v3.00.06 (dot-dashed red line), GiBUU 2019 (dotdot-dashed cyan line), NEUT v5.4.0 (dot-dot-dot-dashed green line), and NuWro 2019 (solid thin purple line).Top: data are shown with total uncertainties, predictions are taken directly from the generators.Bottom: data are shown with shape-only uncertainties, predictions are area-normalized to the data.

Figure 10 .
Figure 10.Single differential cross section as a function of Q 2 (left) and cross section as a function of Eν (right), divided by the GENIE v2.12.2 -NOvA TUNE prediction.Quantities are calculated over the reported muon kinematics space specified by the double-differential measurement.The data are presented statistical and total uncertainties, and compared to the ratio obtained from GENIE v3.00.06 (dot-dashed red line), GiBUU 2019 (dot-dot-dashed cyan line), NEUT v5.4.0 (dotdot-dot-dashed green line), and NuWro 2019 (solid thin purple line).In the top plots, the predictions are unmodified and the error bars on the data represent the total uncertainties, which includes normalization uncertainties that are 100% correlated across bins.In the bottom plots, the predictions are first area-normalized to the data, and the error bars on the data represent shape-only uncertainties, where 100%-correlated normalization uncertainties have been removed.

Table I .
Mass contributions from various elements in the fiducial volume used in this analysis.

Table II .
Fractional uncertainties and correlations across all bins, broken down by source.Averages are taken across all bins reported in this measurement, weighted by the measured cross section, as described in Eqns.3 and 4.

Table III .
Summary of the neutrino event generators and choice of models used by each to generate the inclusive cross-section predictions against which comparisons to this measurement are made.RFG = Relativistic Fermi Gas, LFG = Local Fermi Gas, L-S = Llewellyn Smith, RPA = Random Phase Approximation, R-S = Rein-Sehgal, B-S = Berger-Sehgal, B-Y = Bodek-Yang, PY = PYTHIA, BUU = Boltzmann-Uehling-Uhlenbeck.

Table V :
Double-differential cross-section d 2 σ d cos θµdTµ results table including total and systematic errors cm 2 GeV nucleon × 10 −39 .cos θµ range Tµ range (GeV) Cross section Total Error Stat.Error

Table V :
Double-differential cross-section d 2 σ d cos θµdTµ results table including total and systematic errors cm 2 GeV nucleon × 10 −39 .cos θµ range Tµ range (GeV) Cross section Total Error Stat.Error

Table V :
Double-differential cross-section d 2 σ d cos θµdTµ results table including total and systematic errors cm 2 GeV nucleon × 10 −39 .cos θµ range Tµ range (GeV) Cross section Total Error Stat.Error