Adjusting Neutrino Interaction Models and Evaluating Uncertainties using NOvA Near Detector Data

The two-detector design of the NOvA neutrino oscillation experiment, in which two functionally identical detectors are exposed to an intense neutrino beam, aids in canceling leading order effects of cross-section uncertainties. However, limited knowledge of neutrino interaction cross sections still gives rise to some of the largest systematic uncertainties in current oscillation measurements. We show contemporary models of neutrino interactions to be discrepant with data from NOvA, consistent with discrepancies seen in other experiments. Adjustments to neutrino interaction models in GENIE that improve agreement with our data are presented. We also describe systematic uncertainties on these models, including uncertainties on multi-nucleon interactions from a newly developed procedure using NOvA near detector data.


I. INTRODUCTION
The non-zero value of the reactor mixing angle θ 13 [1][2][3][4] has enabled searches for leptonic CP violation and measurements of the neutrino mass ordering using long-baseline neutrino oscillation experiments with pion-decay-in-flight beams [5][6][7]. Such experiments can also constrain or measure other standard neutrino oscillation model parameters, such as ∆m 2 32 and θ 23 . Long-baseline experiments generally utilize a two-detector design. A smaller near detector (ND) close to the neutrino production target constrains neutrino flux and interaction cross sections. A larger far detector (FD) is positioned to observe the neutrinos after oscillations. Measurements are based on reconstructed neutrino energy spectra observed in the FD, which are compared to simulated predictions for various oscillation parameter values with systematic uncertainties taken into account. ND data are used to adjust FD predictions and constrain systematic uncertainties, via either a simultaneous fit of ND and FD simulation to the respective data samples [8], or by using differences between ND data and simulation to adjust FD simulation [8,9]. In either case, this process relies on simulation to account for oscillations and the differing beam flux and geometric acceptances between the detectors, making the ND constraint on the FD model-dependent. Interactions of neutrinos with nuclei at neutrino energies around 1 GeV, and the resulting final states, are challenging both to describe theoretically and to measure experimentally. As a consequence, systematic uncertainties in neutrino interaction cross sections are typically among the largest uncertainties affecting long-baseline neutrino oscillation measurements, even with the two-detector approach [5,6].
NOvA is a long-baseline neutrino oscillation experiment, utilizing a 14 kton FD located 810 km downstream of the beam source and a functionally identical 0.3 kton ND located approximately 1 km from beam target. The detectors are made from PVC cells of cross section 3.9 × 6.6 cm 2 and of length 3.9 m (ND) or 15.5 m (FD), which are filled with organic liquid scintillator. This results in detectors with 63% active material by mass and a radiation length of 38 cm. Cells are extruded together in units and joined edgewise along the long dimension to produce square planes, which are then stacked perpendicular to the beam direction in alternating vertical and horizontal cell orientations to permit three-dimensional event reconstruction. The near detector additionally has at its downstream end a "muon catcher" composed of a stack of ten sets of planes in which a pair of one vertically oriented and one horizontally oriented scintillator plane is interleaved with one 10 cm-thick plane of steel. Including the muon catcher, the ND can stop muons up to about 3 GeV. The FD is approximately four times longer, wider, and taller than the ND.
High-purity neutrino or antineutrino beams are produced by the NuMI facility at Fermilab [10] according to the current polarity of two magnetic horns that focus and charge-select the parent hadrons. The detectors are located 14.6 mrad off-axis which results in an incoming neutrino energy spectrum narrowly peaked at 2 GeV. This neutrino energy is chosen to optimize sensitivity to oscillations, since ν e appearance and ν µ disappearance probabilities both experience local maxima at an L/E of around 500 km/GeV. The full NOvA experimental setup, including estimates for the neutrino flux from NuMI, is described in Refs. [5,[11][12][13][14][15].
This paper details adjustments made to the neutrino interaction models used in NOvA's simulation and the construction of associated systematic uncertainties. NOvA's 2019 measurements of oscillation parameters [5] use this work 1 .
The data samples and observables used in the analysis are described in Sec. II. Details of the models in the simulation are given in Sec. III and the adjustment procedure is developed in Sec. IV. We discuss systematic uncertainties associated with the adjustments and how we treat them in Sec. V. Finally, we compare our findings to those of other experiments in Sec. VI.

II. DATA SAMPLE AND RECONSTRUCTION
The NOvA data presented here are from a near detector exposure of 8.03 × 10 20 protons on target with the neutrino beam and 3.10 × 10 20 protons on target with the antineutrino beam, totaling 1.48 × 10 6 selected neutrino interactions and 3.33 × 10 5 selected antineutrino interactions. The events used here are the same events selected in the 2019 NOvA ν µ disappearance oscillation results [5], where the selection criteria, efficiencies, and purities are detailed. After selection in the ND, we expect the neutrino beam candidate sample to be composed of 97.1% muon neutrinos and 2.9% muon antineutrinos, with negligible contributions from neutral-current (NC) or other charged-current (CC) neutrino flavors. For the antineutrino beam we expect 90.2% muon antineutrinos and 9.8% muon neutrinos 2 .
Throughout this paper we compare various observables in our data and quantities we derive from them to the predictions we obtain from simulation. For simulated observables, we distinguish the "true," or generated, value from the "reco" value reconstructed in the detector. The energy of muons that stop in the detectors (E µ ) is measured with a resolution of about 3% using track length, while the energy of all other particles, which collectively make up the hadronic recoil system, is measured using calorimetry. For muon neutrino charged-current interactions in NOvA, the visible hadronic energy (E vis had or visible E had ) is the sum of the calibrated observed hadronic energy deposits in scintillator. This is distinct from the fully reconstructed hadronic energy, E had , which also accounts for unseen energy, such as that lost to dead material in the detector or to escaping invisible neutral particles. We measure E had with an energy resolution of about 30%. The reconstructed neutrino energy E ν is the sum of E µ and E had .
The variables E ν , E µ , p µ (the muon momentum), and cos θ µ (the opening angle between the muon and the neutrino beam directions) are estimated as in the ν µ disappearance analysis [15], as is the method for calculating visible E had . We use these, along with the muon mass m µ , to estimate the square of the four-momentum transferred from the initial neutrino to the nuclear system as In conjunction with the energy transfer q 0 , which we measure as E had , Q 2 reco can then be used to approximate the three-momentum transfer as Finally, we combine E had , Q 2 reco , and the proton mass m p to estimate the invariant mass of the hadronic system as

III. SIMULATION
The NuMI beamline, including the 120 GeV protons and the hadrons produced by their interactions, is simulated using Geant4 [16], as is the flux of resultant neutrinos. This neutrino flux is corrected using tools developed by the MINERvA collaboration for the NuMI beam, adding constraints on the hadron spectrum [17]. GENIE 2.12.2 [18,19] (hereinafter referred to as GENIE) is used to predict the interactions of the neutrinos with the detector. Adjustments to GENIE as used by NOvA are the focus of this paper. GENIE prior to our adjustments will be referred to as the "default" simulation.
The simulation of neutrino interactions is separated into distinct parts within GENIE: the initial nuclear state, the hard scatter, and reinteractions of the resultant particles within the nuclear medium. The initial state in the default GENIE configuration is a global relativistic Fermi gas (RFG) model based on the work of Smith and Moniz [20] and modified by adding a high-momentum tail [21] to account for potential shortrange nuclear correlations [22]. GENIE classifies the hard scatter into four primary interaction types. At neutrino energies around 2 GeV, the three most common are: quasi-elastic (QE) interactions, predicted according to the formalism attributed to Llewellyn Smith [23], which result in a single baryon; resonance (RES) processes, predicted according to the model by Rein and Sehgal [24], which result in baryons and mesons via an intermediate hadronic excited state; and what GENIE calls deep inelastic scattering (DIS), predicted using the Bodek-Yang scaling formalism [25] together with a custom hadronization model [26] and PYTHIA6 [27], which results in a broad spectrum of hadrons from inelastic scattering over a large range of hadronic invariant masses. The fourth primary process is the rare instance where a neutrino scatters from the entire nucleus as a coherent whole (COH). As Fig. 1 shows, the default GENIE configuration does not reproduce the visible hadronic energy distribution in the ND neutrino or antineutrino data, undershooting by as much as 25% in the range from 50 to 250 MeV. GENIE, however, does have optional support for the simulation of meson exchange currents (MEC), a process modeled as a neutrino interacting on a nucleon coupled to another nucleon via a meson. Such a process knocks out multiple nucleons from the nuclear ground state in an otherwise QE-like interaction. Two MEC models were available in GENIE that we considered for use, including "Empirical MEC" [28], and the model by the València group (Nieves et al.) [29]. Other models exist but are not implemented in GENIE 2.12.2 [30][31][32]. None of these models explicitly predict the kinematics of the resulting hadrons. Instead, a separate model is necessary to specify how the momentum transfer should be assigned to individual nucleons. The model GENIE uses for all MEC simulation is a so-called "nucleon cluster" model, in which an intermediate nucleon pair whose initial momenta are drawn from the Fermi sea is assigned the total momentum transfer and then allowed to decay isotropically [28].
GENIE also considers final state interactions (FSI) that can occur as the resultant particles traverse the nuclear medium. These are modeled with the hA-INTRANUKE effective cascade model [33,34]. More discussion and further references regarding neutrino-nucleus scattering theory, experiments, and implementation of neutrino interaction software can be found in Ref. [35].

IV. CROSS-SECTION MODEL ADJUSTMENT METHODOLOGY
As each of the interaction types produced by GENIE has independent degrees of freedom and separate uncertainties, it is essential to consider carefully how each model might be adjusted in order to improve data-MC agreement. We first modify the GENIE predictions by incorporating new advances motivated by theory or external data and corroborate them with NOvA ND data in regions where the various modes are expected to be well separated (see Secs. IV A and IV B). After these adjustments, the prediction still disagrees with the ND data, which we attribute to the considerable uncertainty on the spectral shape of MEC events. We reshape and rescale the MEC component so that its sum with the otherwise adjusted simulation matches NOvA ND data, as described in Sec. IV C.
While this procedure explicitly accounts for two-body knockout via MEC, interactions on nuclear pairs formed by short-range correlations between nucleons in the nuclear ground state can also result in a similar final state. The default simulation does not model this explicitly, but our work reshapes the MEC kinematics to match data, effectively adding such missing processes. We thus use the more inclusive term "2p2h" (twoparticle two-hole, describing the ejected particles and the final-state nucleus) to refer to that channel after our model adjustment. 2. Nucleon momentum distribution and long-range nuclear mean field effects in CCQE The more sophisticated local Fermi gas model of the nuclear ground state employed by Nieves et al. [37] predicts a different initial nucleon momentum distribution than the RFG model. This difference, when combined with the effect of Pauli suppression, changes the available kinematic space in QE reactions. Long-range internucleon interactions analogous to charge screening in electromagnetism also modify the kinematics of QE reactions. A popular approach to account for the latter dynamic in calculations uses the random phase approximation (RPA) [37,38]. The combination of these effects significantly suppresses QE reactions at low invariant four-momentum transferred to the nucleus (Q 2 ), and mildly enhances them at higher Q 2 , relative to the RFG prediction. To approximate the result of these two phenomena, we employ the weighting functions based on the València model constructed by MINERvA [39], hereinafter referred to as "QE nuclear model weights." These weights are parameterized in a two-dimensional space of energy and momentum transfer to the nucleus (q 0 , | q|), and are calculated separately for neutrinos and antineutrinos 3 .

Soft non-resonant single pion production
We also reweight GENIE single pion DIS events with invariant hadronic mass W < 1.7 GeV/c 2 to reduce their rate by 57% 4 according to the results of recent reanalysis of bubble chamber data [40]. This is compatible with MINERvA's recent findings [41].
Since that analysis applies only to neutrinos and the analogous GENIE prediction for antineutrinos is very different, we do not apply this correction to antineutrino soft non-resonant single pion production 5 . Similarly, no correction is made to NC channels, as the bubble chamber analysis was for CC channels only.

B. Low-Q 2 resonance suppression
Measurements of neutrino-induced ∆(1232) resonance production [42][43][44][45][46] suggest a suppression at low Q 2 relative to the Rein-Sehgal free-nucleon prediction. Our own ND data reproduces this phenomenon, as seen in the top of Fig. 2. To our knowledge, there is no phenomenology predicting such an effect, though it superficially resembles the effect that the QE nuclear model weights have on QE interactions. We find that applying an alternate parameterization 6 of the QE nuclear model weights discussed in Sec. IV A to RES interactions significantly reduces the tension we observe with our data, as shown in the bottom of Fig. 2.
We therefore reweight all RES events according to this prescription. Formally, the RPA phenomenology may not apply directly to baryon resonance excitation, which requires significant three-momentum transfer to the nucleus even at Q 2 = 0, and thus places RES interactions out of the regime where current RPA calculations apply. However, we employ this procedure as a placeholder for whatever the true effect may be, and invite further input from the theoretical community as to what ingredients may be missing from the model. Data are shown with statistical error bars, while simulation is shown as histograms stacked by interaction type. All cross section adjustments described in this paper are applied, including the addition of the fitted 2p2h described in Sec. IV C, except that the RPA-like low-Q 2 suppression is not applied to RES interactions in the top plots. Neutrino beam is shown at left, antineutrino beam at right.

C. Multi-nucleon knockout (2p2h)
Significant disagreement with the ND data remains even after combining any of the MEC models available in GENIE with the prediction after the modifications described above, as can be seen in Fig. 3. Both the Empirical and Valencia MEC models produce too low of an overall rate, especially at low values of hadronic energy. The visible hadronic energy shapes of Empirical and Valencia MEC are quite different for neutrinos but similar for antineutrinos. It is clear that any MEC model GENIE offers will require significant tuning to reproduce our data. We choose to use the Empirical MEC model as a starting point, as it is the only model available in GENIE that includes a neutral-current component.
The Empirical MEC model is reshaped to create an ad-hoc model that matches data by modifying it in a two-dimensional space of (q 0 , | q|). Simulated GENIE Empirical MEC interactions are divided into 16 bins of energy transfer (from 0 to 0.8 GeV) and 20 momentum transfer bins (from 0 to 1 GeV/c). Of these, 120 bins are kinematically disallowed. Scale factors for each of the remaining 200 bins in (q 0 , | q|) are incorporated as Gaussian penalty terms into a χ 2 fit, each with 100% uncertainty. For this fit, the non-2p2h portion of the simulation is adjusted as described in this paper, and the 2p2h component is reweighted as dictated by the penalty terms. A migration matrix is used to convert the (q 0 , | q|) prediction into a binned 20x20 space of visible hadronic energy E vis had (from 0 to 0.4 GeV) and reconstructed three-momentum transfer | q| reco (from 0 to 1 GeV/c). This prediction in reconstructed variables is then compared to the ND data in the fit. The small (2%) antineutrino MC component is left in its default state when fitting the neutrino beam simulation to data. The process is repeated for the antineutrino beam data and MC, except in this case the 2p2h fit for neutrinos is applied first to the larger (about 10%) neutrino component in the antineutrino beam MC.
The resulting weights are shown in Fig. 4. Since true q 0 and E vis had are strongly correlated variables, the enhancement of events at low values of q 0 compensates for the deficit of simulated events at low visible hadronic energy seen in Fig. 3. In the antineutrino beam sample there is less discrepancy at low E vis had than in the neutrino beam sample, and thus the antineutrino weights show a smaller enhancement at low q 0 . Additionally, events in the higher q 0 tail are suppressed for antineutrinos. These features are evident in Fig. 5, which compares the unaltered Empirical MEC distributions in energy transfer and momentum transfer to the reweighted distributions.   Default"), the default configuration with the addition of unadjusted Empirical MEC ("+MEC"), and after all the adjustments of Sec. IV ("Final"). The "Before selection" column indicates the fully adjusted fractions before selection, illustrating the important role acceptance and reconstruction efficiencies play in the ND. Fractions may not add to precisely 1.00 due to rounding.

D. Summary of adjustments to central value prediction
In summary, the NOvA prescription for adjusting GENIE cross-section models to incorporate external data constraints and to improve agreement with NOvA ND data is to start with GENIE, using the Empirical MEC model, and: 1. Change CCQE M A from 0.99 to 1.04 GeV/c 2 ; 2. Apply València nuclear model weights from MINERvA, using the (q 0 , | q|) parameterization for QE and the Q 2 parameterization for RES; 3. Apply a 57% reduction to soft non-resonant single pion production events from neutrinos; 4. Apply separate ν andν weights in (q 0 , | q|) derived from NOvA ND data to Empirical MEC interactions.
The effect of each step is shown in Fig. 6. The default GENIE simulation has a large deficit of events in the MEC region in both beams when compared to data, though the neutrino beam prediction has a 5% excess in the lowest hadronic energy bin. The QE modifications particularly affect the low E vis had region due to the suppression from the nuclear model. The adjustment to RES and DIS widens the deficit, then by design the 2p2h fit modifies the shape of this component to improve agreement. The predicted composition of the sample before and after the tuning procedure is given in Table I.
The final distributions of E vis had and | q| reco after all adjustments are shown in Fig. 7. The modified simulation largely matches data (by construction) in regions where 2p2h is significant. The lowest visible hadronic energy bin in both beams still shows disagreement, mostly due to smearing from the quantities being modified (q 0 , | q|) to the reconstructed quantities (E vis had , | q| reco ) used in the fit. There are residual discrepancies in the regions dominated by pion production, which suggests further model adjustments may be warranted. Figure 7 also shows the final neutrino energy distribution, which is the key variable in neutrino oscillation measurements. The shape of this distribution, and the resolution with which NOvA measures it, is largely unchanged by the adjustment procedure, since the NOvA detectors are calorimeters and the changes do not significantly change the amount of invisible energy. According to the simulation, the mean bias (E reco ν − E true ν )/E true ν is -3.6% (-2.5%) for neutrinos (antineutrinos) with GENIE's default prediction and -2.3% (-2.1%) after all the adjustments; the RMS of this variable shifts from 10.6% (9.3%) to 10.5% (9.3%). 7 Figure 8 shows the visible hadronic energy in bins of momentum transfer, illustrating that the adjusted 2p2h component resides at intermediate values of q 0 and | q|, as expected from observations in electron scattering [47] and in MINERvA [48,49]. This is a key indicator that the discrepancy between the default simulation and ND data is likely due largely to 2p2h interactions. Other kinematic distributions comparing data and simulation can be found in Appendix A.

V. CROSS-SECTION SYSTEMATIC UNCERTAINTIES
GENIE includes an evaluation of many cross-section uncertainties and enables corresponding adjustments to model parameters. We employ this uncertainty model, the details of which can be found in the GENIE manual [19], largely unchanged. However, we substitute our own treatment in several instances where different uncertainties are warranted, as described in the following sections.

A. Quasi-elastic interactions
The default GENIE systematic uncertainty for CCQE M A is +25%/-15%. This uncertainty was constructed to address the historical tension between bubble chamber and NOMAD measurements [50], and MiniBooNE [51], tension which is now largely attributed to be due to multi-nucleon effects [52]. As we explicitly add these multi-nucleon effects and their associated uncertainties separately, we reduce the size of the CCQE M A uncertainty to 5%, which is a rough estimate of the free nucleon scattering uncertainty derived from bubble chamber measurements [53][54][55][56][57].
In addition to the central value weights discussed in Sec. IV A.2, the València CCQE nuclear model weights supplied by MINERvA include separate sets of weights that (when applied to the GENIE RFG distributions) produce alternate predictions for the València model under enhancement and suppression uncertainties [39]. Separate weights are generated for neutrinos and antineutrinos. We include these variations in the uncertainties we consider.

B. Resonance interactions
As discussed in Sec. IV B, the Q 2 parameterization of the QE nuclear model effect applied to RES is a placeholder for an unknown effect. Therefore, we take a conservative 100% one-sided uncertainty on this correction. This permits the effect to be turned off, but not increased, and it cannot change sign. This is the largest systematic uncertainty in NOvA's measurement of θ 23 [5], and is correlated between neutrinos and antineutrinos.

C. Deep inelastic scattering
GENIE's uncertainty model includes uncorrelated 50% normalization uncertainties for DIS events with one-or two-pion final states (any combination of charged or neutral) and W < 1.7 GeV/c 2 . 8 However, there is no corresponding normalization uncertainty for DIS with W > 1.7 GeV/c 2 , or for any events with pion multiplicity larger than two. Moreover, the sharp discontinuity going from 50% to 0% when crossing the W = 1.7 GeV/c 2 boundary leads to unphysical variations when used to produce alternate predictions. We therefore replace the low-W GENIE DIS normalization uncertainties with 32 [4 (0π, 1π, 2π, > 2π) × 2 (CC, NC) × 2 (neutrino, antineutrino) × 2 (interaction on neutron, proton)] independent, uncorrelated 50% normalization uncertainties for all DIS events up to 3 GeV/c 2 in W . These uncertainties drop to 10% for the W > 3 GeV/c 2 region, which is more consistent with previous measurements of DIS at higher energy 9 . A comprehensive summary of the available data and corresponding theory is given in Ref. [58].

D. 2p2h
We include three types of 2p2h uncertainty, all of which we take as uncorrelated between neutrinos and antineutrinos, for a total of 6 independent uncertainties. Throughout, we neglect the influence of shortrange correlations on the uncertainties we consider since the 2p2h contribution to the neutrino interactions considered in this work is expected to be dominated by MEC [59].

Target nucleon pair identities
A CC MEC interaction always involves a target nucleon whose identity (proton or neutron) is dictated by charge conservation. The identity of the second nucleon, coupled to the first in the interaction, is determined by the model. We examine various theoretical models to determine the relative proportions of neutrons versus protons in the struck (initial state) nucleon pairs and use these predictions to construct an uncertainty. For neutrinos, we are interested in the fraction of target pairs that are neutron-proton, R ν = np/(np + nn), which for the València model included in GENIE averages 0.67 over the kinematic range of interest. A detailed study during the development of the SuSA MEC model concluded that, over a range of kinematics, their fraction is 0.8-0.9 [31]. The Empirical MEC model in GENIE defaults to a value of 0.8. Though the València model predicts R as a function of the momentum transfer, Empirical MEC does not, and we do not have a full simulation of the SuSA model to study the impact in our phase space. For this analysis we therefore retain R ν = 0.8 as a fixed central value and take the range 0.7-0.9 as a 1σ uncertainty. In future work we plan to study the effect of the differing models' predictions as a function of kinematics in more detail. For antineutrinos, we use the same central value and uncertainty range for the Rν = np/(np + pp) ratio. This uncertainty has a small effect on predictions of observables; the expected visible hadronic energy shapes of R = 1 vs R = 0 events are shown in Fig. 9.

Energy dependence of total cross section
The second uncertainty addresses the difference between MEC models in cross section as a function of neutrino energy. Four MEC models are examined: Empirical [28], València [29], that of the Lyon group (Martini and Ericson) [30], and SuSA [31]. As our tuning procedure enforces a normalization inferred from our data, we are concerned mostly with shape differences; therefore, we rescale the predictions. In principle, we prefer to normalize at higher energies where the predicted spectra flatten, but several models do not extend this far. Thus, we take the following approach: the València prediction from GENIE is scaled to match Empirical MEC at 10 GeV; the SuSA prediction is scaled to match Empirical MEC at 4 GeV (the highest-energy prediction in [31]); and the Lyon prediction is scaled so that its peak is the same as that of Empirical MEC. Our rescaled predictions for σ(E) from the models are shown in Fig. 10a. We compute the ratios of the renormalized model predictions to Empirical MEC and construct a function which approximately envelopes the variations, as shown in Fig. 10b. This function becomes an energy-dependent 2p2h normalization uncertainty.
This procedure is based on neutrino MEC models. Since fewer models that consider antineutrinos are available, the same envelope is used (uncorrelated) for antineutrinos. 3. 2p2h dependence on non-2p2h prediction The 2p2h fit reshapes the Empirical MEC interactions such that the total simulation will match ND data. Any imperfections in other parts of the simulation will consequently be absorbed into the resulting 2p2h sample. We can quantify this uncertainty by examining the dependence of the 2p2h fit on other systematic uncertainties. These reactions are known to occupy a region of energy transfer in between QE interactions (at low q 0 ) and RES interactions (at higher q 0 ); this holds true in E vis had as well. In general, uncertainties that affect the E vis had distribution of the non-2p2h prediction shift the mean to be higher or lower in q 0 , and thus more like a purely RES or QE spectrum. As a result, the fitted 2p2h spectrum moves in the opposite direction in q 0 . A similar effect holds in | q|. Using the largest non-2p2h cross-section systematic uncertainties, we apply correlated 1σ shifts to create the largest q 0 -shifting distortions allowed by our uncertainty treatment, which conservatively bound this effect.
The shifts listed in Table II are combined to distort the non-2p2h simulation to be more more "RES-like" or "QE-like", resulting in a fitted 2p2h prediction that is more "QE-like" or "RES-like" respectively. The uncertainties in the table are either standard GENIE systematic uncertainties or are described herein.
TABLE II. Correlated systematic uncertainty shifts used to make the non-2p2h simulation more "RES-like" or "QElike" before fitting the 2p2h component.
The 2p2h fitting procedure is carried out in each of these two scenarios, for both neutrinos and antineutrinos separately, to create ±1σ shape uncertainties. The differences in the fitted q 0 predictions are illustrated in Fig. 11. We anticipate that 2p2h predictions made using these alternative underlying model assumptions should bracket the unknown true 2p2h response. In the future we anticipate considering other potential sources of 2p2h uncertainty that we have assumed to be subdominant here, including the assignment of final-state energies to the nucleons in the nucleon cluster model in GENIE.

E. Summary of cross-section model uncertainties
Our modifications and additions to the default GENIE model uncertainties are summarized below. In this section, "uncorrelated" means that parameters in the uncertainty are allowed to vary separately for neutrinos and antineutrinos; "correlated" indicates that neutrinos and antineutrinos use the same values.
We also introduce three additional uncertainties: 1. QE nuclear model uncertainties (different for neutrino and antineutrino; uncorrelated); 2. A 100% uncertainty on the RES low-Q 2 suppression, which can never go above 100% or negative (correlated); 3. Three 2p2h uncertainties: one covering uncertainty in target nucleons, one addressing uncertainties in the energy dependence of the cross section, and one treating uncertainties in the (q 0 , | q|) response (all uncorrelated).
The combined cross-section uncertainties are shown in Fig. 12. The adjusted neutrino simulation agrees with data somewhat better than the antineutrino simulation, but in both cases the data lies within the uncertainty band.

VI. COMPARISONS TO OTHER OBSERVATIONS
As shown in Fig. 7 and Appendix A, the total inclusive prediction, including the 2p2h component tuned in (q 0 , | q|) space and fit in (E vis had , | q| reco ), can reproduce our observed ND distributions in numerous kinematic variables. MINERvA, an on-axis experiment using the same neutrino beam as NOvA, has performed an analogous 2p2h tuning procedure with their inclusive neutrino-mode data set [48]. They use GENIE with the same QE nuclear model weights described in Sec. 2, and apply a correction to non-resonant single pion production similar to that in the NOvA prescription, but use the València MEC model. In their procedure, the values of a two-dimensional Gaussian are taken as weights to the MEC prediction, and the Gaussian's parameters are fitted in order to match the observed distributions [60]. They find good agreement with their antineutrino data using this adjusted model with no further modifications [49]. The result of replacing the 2p2h component of the NOvA fully adjusted simulation with the MINERvA tuned 2p2h prediction is shown in Fig. 13. Qualitatively, the MINERvA model results in a similar overall prediction to the NOvA model, mostly falling within the 1-σ uncertainties.
The T2K collaboration uses NEUT [61,62] instead of GENIE to simulate neutrino interactions for their primary neutrino oscillation analysis. In their recent work [6] they also use implementations of the València models for the central value prediction of both QE and MEC processes. Among the uncertainties they consider for QE is a parameterized version of the nuclear model calculations for long-range correlations that is similar to that used by NOvA and MINERvA. Uncertainties in the MEC model are bounded between two extreme cases: a prediction using only those MEC diagrams coupling to a ∆-resonance, and a prediction removing all the ∆ channels. The T2K fit pulls this 2p2h shape uncertainty to the maximum allowed value [63]. The 2p2h normalization is also pulled to be 50% larger than the default prediction. This is consistent with the findings by NOvA and MINERvA that using an unaltered version of the València model is insufficient to describe data.

VII. CONCLUSIONS
We find that modifications to the default GENIE 2.12.2 model significantly enhance the agreement between selected muon neutrino candidates in the NOvA ND data sample and simulation across a variety of kinematic variables. Corrections to the QE and soft non-resonant single pion production predictions based on reevaluated bubble chamber measurements are included. Improved nuclear models are also used to adjust the kinematics of QE scattering. Furthermore, suppression at low Q 2 on resonant pion production is imposed as supported by observations in other experiments and our own ND data. The Empirical MEC model in GENIE is tuned to match data in our ND. A set of systematic uncertainties are created, addressing potential weaknesses in the models and bounding the results of our own tuning procedure with ND data.
We will continue to incorporate constraints from other measurements as well as advances in cross-section modeling into our predictions and reduce the impact of systematic uncertainty on our analyses. Such improvements will not only benefit NOvA and other current experiments, but will be necessary for future experiments such as DUNE, which has stringent requirements on its systematic uncertainty budget [64]. In each bin, the 1σ deviations from nominal for each cross-section uncertainty are added in quadrature to obtain the band, which has significant bin-to-bin correlations.  FIG. 13. Comparison of reconstructed visible hadronic energy distribution in ND data (black dots) to various simulations for neutrino beam (left) and antineutrino beam (right) running. The solid black curves correspond to GENIE predictions with the full set of adjustments described in this paper, while the red and purple dotted curves are the simulation with +1 and −1σ shifts from the 2p2h (q0, | q|) response systematic uncertainty shown in Fig. 11, respectively. Also shown in solid blue is the result of replacing our tuned 2p2h with MINERvA's tuned 2p2h prediction. The shaded gray histogram represents the GENIE prediction for non-2p2h interaction channels.