Parametrized uncertainties in the spectral function model of neutrino charged-current quasielastic interactions for oscillation analyses

A substantial fraction of systematic uncertainties in neutrino oscillation experiments stem from the lack of precision in modeling the nuclear target in neutrino-nucleus interactions. Whilst this has driven significant progress in the development of improved nuclear models for neutrino scattering, it is crucial that the models used in neutrino data analyses be accompanied by parameters and associated uncertainties that allow the coverage of plausible nuclear physics. Based on constraints from electron scattering data, we develop such a set of parameters, which can be applied to nuclear shell models, and test their application to the Benhar et al spectral function model. The parametrization is validated through a series of maximum likelihood fits to cross-section measurements made by the T2K and MINERvA experiments, which also permit an exploration of the power of near-detector data to provide constraints on the parameters in neutrino oscillation analyses.

A substantial fraction of systematic uncertainties in neutrino oscillation experiments stem from the lack of precision in modeling the nuclear target in neutrino-nucleus interactions.Whilst this has driven significant progress in the development of improved nuclear models for neutrino scattering, it is crucial that the models used in neutrino data analyses be accompanied by parameters and associated uncertainties that allow the coverage of plausible nuclear physics.Based on constraints from electron scattering data, we develop such a set of parameters, which can be applied to nuclear shell models, and test their application to the Benhar et al. spectral function model.The parametrization is validated through a series of maximum likelihood fits to cross-section measurements made by the T2K and MINERvA experiments, which also permit an exploration of the power of near-detector data to provide constraints on the parameters in neutrino oscillation analyses.

I. INTRODUCTION
Modeling neutrino-nucleus interactions remains a major challenge for current and future neutrino oscillation experiments [1].Whilst statistical uncertainties remain large in the ongoing accelerator-based long-baseline neutrino oscillation experiments, T2K [2] and NOνA [3], the future experiments, Hyper-Kamiokande (HK) [4] and DUNE [5], are likely to be dominated by uncertainties related to neutrino interactions.To tackle this issue, experiments are designing and constructing more sophisticated near detectors [4,6,7].However, more accurate neutrino interaction models are a critical ingredient to precisely extrapolating new near-detector measurements into robust constraints on the uncertainties in neutrino oscillation analyses.Beyond just incorporating better models, it is also necessary for future neutrino oscillation analyses to include a comprehensive set of uncertainties to cover their plausible variations.Establishing such lists of uncertainties remains one of the most significant challenges for robustly reducing systematic uncertainties using near-detector data in neutrino oscillation analyses.
At E ν ∼ 1 GeV energy, neutrino interactions with nucleons bound within nuclear targets are significantly impacted by nuclear effects which can have a substantial impact on the interaction cross section and on the kinematics of the final-state particles.These can affect the bias in metrics for reconstructing neutrino energy from observable final-state interaction products which, if modeled incorrectly, directly biases the measurement of neutrino oscillation parameters [1].In particular, due to the removal energy and the Fermi motion of the bound nucleon, the neutrino interacts on an off-shell, non-static particle which can significantly alter interaction kinematics with respect to the case of a free-nucleon target.Multiple processes can further affect the final-state kinematics: this is the case for Pauli blocking (PB), which prevents some interactions to occur, and for final-state interactions (FSI), which describe how the interaction between the outgoing hadrons and the residual nucleus affects the process.
A relatively sophisticated way to model chargedcurrent quasielastic interactions (CCQE)-the dominant interaction channel for T2K and HK-is provided by the spectral function model (SF) by Benhar et al. [8].It is a nuclear shell model built largely from electron scattering data, and allows a detailed description of the initial nuclear state.Through its implementation in the NEUT neutrino-nucleus interaction event generator [9], the SF model is used in the T2K experiment's recent measurements of neutrino oscillations [10] and neutrino-nucleus cross sections [11][12][13][14].It has also been used to study the sensitivity of T2K's near-detector upgrade [15].
Shell-based models lend themselves to consistently defining an ensemble of uncertainties which affect outgoing lepton and hadron kinematics in neutrino interactions.The shell structure presents natural degrees of freedom that can be varied, and their determination from electron scattering analyses can motivate possible model variations.In this paper, we detail a set of parameters to alter the SF alongside associated uncertainties.The general scheme is applicable to any shell-based model.The parameters applied to the SF model are benchmarked by fitting them to neutrino-nucleus cross-section measurements from the T2K and MINERvA experiments.
The paper is organized as follows: in Sec.II the Benhar et al.SF model is summarized, and Sec.III introduces a parametrization of systematic uncertainties for the SF model as well as for other physics processes that affect comparisons with cross-section measurements.In Sec.IV we show how the use of this parametrization in fits of the NEUT generator to existing neutrino cross-section measurements can improve the agreement of the SF model with recent measurements from the T2K and MINERvA experiments.Means of mitigating effects from "Peelle's pertinent puzzle" (PPP) [16,17] when fitting the measurements are also presented and applied.In Sec.V a qualitative discussion of the parametrization's impact on neutrino oscillation analyses is presented, alongside an assessment of the scope for possible constraints from near detectors.Finally, the discussion and conclusions are presented in Sec.VI.

II. THE BENHAR SPECTRAL FUNCTION MODEL FOR NEUTRINO OSCILLATION EXPERIMENTS
Neutrino experiments rely on models available in Monte Carlo event generators, such as NEUT [9], NuWro [18] and GENIE [19], to provide a model of neutrinonucleus interactions for neutrino oscillation measurements.Various prescriptions describe the initial state of the bound nucleons within a target nucleus.Generally speaking, neutrino oscillation experiments use either a relativistic global Fermi gas (RFG) model, a local Fermi gas (LFG) model, or a spectral function (SF) model.The RFG is a simplistic approach that describes the nuclear ground state and Fermi motion, which has been commonly-used to model GeV-scale neutrino scattering.It considers the target nucleons as independent fermions within a constant binding potential.The LFG model is a more realistic approach that introduces the radial dependence in the nuclear potential under the local density approximation (LDA).The SF model provides a sophisticated description of the nuclear ground state, featuring a shell structure of the nucleus as observed in electron scattering data, with an additional theory-driven component to describe "short-range correlations" (SRC) between nucleons within the nucleus.Fig. 1 shows a comparison of the predictions between the three models for the distribution of the nucleon removal (or missing) energy, E m , and the initial-state nucleon (or missing) momentum, p m , in carbon.Closely following the prescription in Ref. [20], for CCQE neutrino-nucleus interactions with a single nucleon in the final state before consideration of FSI (i.e.ν + A → ℓ + N + A ′ , where A and A ′ are the initial-state nucleus and the final-state remnant respectively and N is the outgoing nucleon), these can be defined as1 : where E ν and ⃗ p ν (E ℓ and ⃗ p ℓ ) are the incoming neutrino (outgoing charged lepton) energy and momentum respectively; T N and p N are the kinetic energy and momentum of the pre-FSI outgoing nucleon respectively; T A ′ is the kinetic energy of the nuclear remnant; and ∆m N accounts for the mass difference of the initial-and finalstate nucleon.T A ′ is calculated from p N assuming a mass of the ground state of the remnant nucleus following the interaction.
Fig. 1 clearly highlights the nuclear shell structure of carbon encoded in the SF model, with a sharp p-shell at E m ∼ 18 MeV and a diffuse s-shell at E m ∼ 35 MeV, which is not captured by the Fermi gas-based models.In Fermi-gas models the kinetic energy from the Fermi motion of the stuck nucleon contributes to overcoming the removal energy which gives the parabolic shapes observed.Both SF and LFG have a minimum removal energy which is derived from electron scattering data or the difference between initial state and remnant nucleus masses respectively (for the LFG case, more details can be found in [21]).
The distribution of the initial-state nucleon momentum is also different between the models, as displayed most clearly in the left panel of Fig. 2. The RFG model assumes that p m is uniformly distributed in 3-momentum, ⃗ p, below the Fermi momentum p F , which yields a cliff feature at p m ∼ 220 MeV/c.The LFG and SF models both predict narrower p m distributions than the RFG model.A high momentum tail appears only in the SF model, a direct consequence of the contributions from SRCs.
These differences in nuclear ground-state modeling can manifest themselves in important ways for neutrino oscillation analyses.Of particular importance is the impact they can have in determining the bias of neutrino energy reconstruction.The oscillation analyses in T2K and HK focus on single-track events, meaning the selections are dominated by CCQE, two-particle two-hole (2p2h), and single-pion interactions where the pion is missed or absorbed.In single-track events, the incoming neutrino's energy is estimated using only the kinematics of the outgoing charged lepton by assuming that the interaction is CCQE on an initial-state nucleon at rest, with an average binding energy E b , and that there is a single outgoing nucleon2 , where m N is the mass of the struck nucleon, mN = m N − E b , and m l , E l , p l , and θ l are the mass, the energy, the momentum, and the angle between the charged lepton and the neutrino, respectively.The right panel of Fig. 2 shows the bias of this estimator with respect to the true neutrino energy, E ν , for the different groundstate models.The smearing is dominated by the isotropic Fermi motion of the bound nucleon.The simplistic RFG model is offset due to its shifted E m distribution with respect to the other models, as shown in Fig. 1, while the LFG and SF models are more similar.Since the SF model accounts for the SRC contribution, E QE ν tends to underestimate the neutrino energy for these events, hence the longer tail in the negative region.
The SF model, which is the focus of this work, provides a description of the two-dimensional distribution of the missing energy, E m , and the missing momentum, p m , of the nucleons within the nucleus.It is implemented in both the NEUT and NuWro event generators, where the joint probability distribution of E m and ⃗ p m of nucleons within the nucleus is written as: where P MF (⃗ p m , E m ) is the mean field (MF) part and P corr (⃗ p m , E m ) corresponds to the contribution from SRC nucleons.The MF term describes the nucleus as a shell model and is tightly constrained with electron scattering data.Within neutrino event generators, SRC come from generally high p m (above the RFG Fermi level) CCQE interactions, modeled as if with a single "target" nucleon but where a "spectator" nucleon is also present.The spectator has similar magnitude but opposite direction ⃗ p m to the target nucleon.Prior to interaction, the spectator and target only remain inside the nucleus because of their correlation maintaining a collectively bound state.
Once the target nucleon is ejected from the nucleus, the spectator, despite not receiving any of the interaction's momentum transfer, has too high p m to remain bound without its pair and so also leaves the nucleus.It should be noted that SRCs are implemented in generators as being distinct from 2p2h, which directly consider the cross section from a series of Feynman diagrams in which the neutrino interaction takes place with two nucleons bound by meson exchange currents and that the momentum transfer of the interaction is shared between them [22,23].Whilst the original SF model calculates the MF and SRC components of the SF separately, the implementation in event generators is as a single integrated SF which the generators then choose to factorize into an MF and SRC parts as disjunct areas in E m , p m space.The two-dimensional (p m , E m ) distributions for oxygen and carbon are shown in Fig. 3, where the MF and SRC regions as defined in NEUT (solid) and NuWro (dashed) are also drawn.
Although SF represents one of the more sophisticated nuclear ground-state models implemented in neutrino event generators, it is still, like most other implemented models, based mostly on the plane-wave impulse approximation (PWIA).PWIA assumes that the interaction cross section can be calculated as the incoherent sum of scattering from free nucleons with a given initial energy and momentum, and that the outgoing nucleon is not affected by the nuclear potential (i.e. it exits as a plane wave).In other words, only one nucleon of the target nucleus is involved in the interaction, and there is no consideration of FSI as the nucleon propagates through the nucleus at the level of the cross-section calculation.
Within the SF model's implementation in neutrino event generators, several effects beyond PWIA are considered using ad-hoc approaches.These include Pauli blocking (PB) and FSI via an intranuclear cascade [24] and/or a correction to the cross section based on the distortion of the outgoing nucleon via a nuclear potential (see e.g.[25,26]).PB accounts for the statistical correlations between the struck nucleon and the remnant nucleus.A simplistic approach to estimate its effect, which is used in NEUT, consists of modifying the spectral function by: where θ is the Heaviside step function, ⃗ q is the transferred momentum, and p F is the Fermi surface momentum, which is nucleus dependent.This amounts to considering that all the states below Fermi momentum are occupied, which suppresses a portion of the available phase space for the outgoing nucleon.Other prescriptions are also possible, such as the LDA-based PB implemented in NuWro.
PB is a first-order effect beyond the impulse approximation part of PWIA.Whilst PWIA is justifiable for high energy transfer, the interactions with low energy transfer (q 0 ≲ 70 MeV)-which constitute a significant fraction (∼ 20%) of CCQE neutrino interactions in oscillation experiments-need to account for other effects beyond the PWIA.In the SF model, one such effect is the aforementioned impact of FSI on the cross section.The impact of FSI on outgoing nucleon kinematics is usually included in event generators via hadron transport models such as intranuclear cascades [24] with re-interaction probabilities derived from data.Another approach to considering the impact of FSI on an SF model is presented in Ref. [27], where modifications to the inclusive cross section are calculated by considering the outgoing nucleon in a nuclear optical potential.The relationship between these two approaches to describing FSI is described in Ref. [28].Within event generators, NEUT only implements an intranuclear cascade whilst NuWro offers both a cascade and an optional optical potential correction (which is not applied in Fig. 3).
Recent neutrino-nucleus cross-section measurements provide a crucial benchmark for all these models.Whilst CCQE-focused measurements which integrate over nucleon kinematics often favour Fermi gas-based models, thanks to the strong "RPA" correction to account for correlations between initial-state nucleons applied as a step away from the impulse approximation [29], measurements that include nucleon kinematics tend to disfavor them and prefer the SF model due to its predictive power of the outgoing nucleon kinematics (although it should be noted that no model is able to achieve reasonable agreement with all data) [30,31].For instance, the singletransverse variables [32] allow a probe of nuclear effects, particularly with the measurement of the transverse momentum imbalance δp T defined from the kinematics of the outgoing charged lepton and nucleon as: where ⃗ p T l and ⃗ p T N are the momenta of the charged lepton and the struck nucleon projected on the transverse plane with respect to the incident neutrino direction, respectively.It corresponds to the transverse projection of the Fermi motion of the initial-state nucleon, modified by nuclear effects such as FSI.Measuring δp T therefore probes the nuclear effects experienced by the struck nucleons.The bulk of its distribution is mainly due to Fermi motion, while its tail is more sensitive to FSI and multinucleon processes.The T2K [11] and MINERvA [33] measurements of δp T show that the SF model is indeed more suitable than Fermi-gas models to describe their measurements (although it remains unable to give a quantitatively good description of all the measurements).Recent analyses from MicroBooNE [34] offer new measurements, which have not yet been confronted with SF models.

III. SYSTEMATIC UNCERTAINTIES
As discussed in Sec.II, the SF model offers a more detailed description of the nuclear medium in comparison with widely-used alternative models like RFG or LFG.For this reason, the T2K experiment opts to use the SF model as a reference for CCQE interactions [10].As described in Sec.I, a detailed investigation of the potential uncertainties related to the SF (or PWIA-based shell models in general) and the development of a comprehensive associated uncertainty parametrization are of great interest.Here, we present a parametrization of systematic uncertainties that are subsequently benchmarked in Sec.IV via fits to cross-section measurements.
In addition to defining an uncertainty parametrization, prior values for the parameters' constraints are also assigned for fits in Sec.IV-from both experimental data and theoretical arguments.In general, we opt for loose conservative priors, thereby allowing the data to dominate the constraints on most of the parameters in the fits.This approach is motivated by the fact that priors from experimental data might not be entirely applicable for neutrino scattering (e.g.taking priors from electron scattering data) or that the prior constraints are based on model-dependent assumptions.The exact priors used are stated and justified in the following sub-sections.Overall, Tab.II summarizes all the parametrized systematic uncertainties discussed in this section, along with their central values and prior uncertainties.
A. Shell-model uncertainties Nuclear shell models provide several natural degrees of freedom that can be parametrized and varied as systematic uncertainties.The SF model contains several discrete shells with a particular size and shape in addition to an SRC component.The occupancy and shape of each shell as well as the relative strength of the SRC component provide such natural freedoms.Moreover, the degree to which these can vary can be inspired by the electron scattering data that informs the nominal SF contribution.

Shell occupancy
As shown in Fig. 3, the missing energy distributions in the SF model exhibit multiple peaks, which correspond to the nucleon energy levels in a shell model.Electron scattering measurements show that the relative strength of each shell may differ from those predicted by a simple independent-particle shell model [25,26].To account for this, a parameter that modifies the shell occupancy is introduced as a normalization uncertainty for each shell in the MF region of the SF model.Weights are assigned to CCQE events in this region following: where N shell is the normalization parameter of a given shell, and E shell and σ shell correspond to the center and the width of the Gaussian function, g shell (E m ), which are fixed for each shell.In total, this gives two shell normalization parameters for interactions with carbon and three for interactions with oxygen.The fixed values of E shell and σ shell are derived from an analysis of the missing energy distributions in NEUT, and are shown for each shell in Tab.I.The impact of varying N shell for the carbon pand s-shells is shown in Fig. 4.
As can be seen in the E m distribution from Fig. 4, this approach results in a thinner or thicker shell for higher or lower N shell values respectively and allows the normalization of the cross section to change (i.e. the total spectral function normalization is not scaled down to account for an increase in a particular N shell ).Note also that there is no limit to the range of E m each shell normalization uncertainty dial is applied to.Each event gets weights from shell occupancy dials for every shell, but if an event is far away from a shell center then the corresponding shell occupancy weight would require very large values of N shell to give the event a significant weighting.One important effect of the shell occupancy parameters is how they alter the total p m distribution, since the initial-state nucleon's momentum distribution differs between the shells, as shown in Fig. 5. Therefore, a change in the relative strength of a shell impacts the shape of the overall distribution of p m , and consequently the distribution of δp T , as is illustrated in the lower panel of Fig. 4.Such freedom allows for variation of the outgoing lepton and nucleon kinematics in a way that is consistently propagated through the model.
As discussed, we opt to keep the prior uncertainties on the N shell parameters conservative.The chosen values are reported in Tab.I.These values cover differences in the shell occupancy of the SF model (derived from (e, e ′ p) data) to those expected in an independent particle shell model [25,26], whilst ensuring the same variation for each parameter changes the total CCQE cross section by the same amount (∼ 10%).Fig. 6 shows the impact of varying each shell normalization parameter N shell on the total CCQE cross section, and Fig. 4 demonstrates potential sensitivity to variations through measurements of δp T shape.If the parameters are pulled far outside the prior uncertainties, this could indicate that the parameters are acting as effective degrees of freedom to account for physics beyond the PWIA.

Missing-momentum shape
Electron scattering data from Ref. [20] provides measurements of the missing momentum distributions p m for carbon.By comparing the predicted NEUT distributions of p m in each shell (which are taken directly from Benhar SF model predictions) with this data, small shape differences appear, as illustrated in Fig. 7.This shows two "extreme" distributions of the measured missing momentum which correspond to (e, e ′ p) kinematics with the most different four-momentum transfer, Q 2 , compared to the NEUT prediction.This builds a missing-momentum shape uncertainty for each shell that changes the shape of the p m distribution.Each shell parameter is defined between -1 and 1, corresponding to the two extreme p m distributions, and a linear extrapolation is implemented beyond the range [−1, 1].These p m shape parameters are expected to mainly affect experimental observables sensitive to the initialstate nuclear momentum.This is demonstrated in Fig. 8, which shows the impact of the p-shell shape parameter on the distributions of the neutrino-muon angle and the transverse momentum imbalance.The lepton kinematics have no discernible sensitivity to these variations, whereas the bulk of δp T is affected by this uncertainty since it corresponds, in the absence of other nuclear effects, to the transverse projection of the Fermi motion.

Short-range correlations
As previously discussed in Sec.II, the SRC contribution implemented in NEUT and NuWro neutrino event generators corresponds to CCQE events yielding two outgoing nucleons following the primary interaction.Calculations performed by Benhar et al. [8] are implemented in neutrino event generators via tables of the total SF distribution containing both the MF and SRC components.Whilst these are separate components of the original SF calculation, the information on their shape is not preserved in the generator implementations.Therefore, it is necessary that generators develop a scheme to separate the two components.NEUT uses hard cuts on E m and p m to distinguish between them as shown in Fig. 3: an SRC two-nucleon knock-out only occurs if the neutrino interacts with a nucleon for which E m > 100 MeV or p m > 300 MeV/c.The spectator nucleon of the pair is taken to have opposite isospin and momentum compared to the "active" interacting nucleon.Using this implementation, NEUT predicts that interactions with SRC pairs represent ∼ 5% of the total CCQE interactions for both carbon and oxygen.Alternatively, the SF implementation in NuWro takes a different approach by making nonrectangular cuts in the (p m , E m ) phase space adapted to each target in a more phenomenological manner, also shown in Fig. 3. Furthermore, while the hard cuts in NEUT fully determine the MF-SRC separation, NuWro applies an additional condition to allow for the knock-out 6.The impact of the shell normalization parameters (N shell ) on the total CCQE cross section when applied to the SF model implemented in NEUT for each shell in carbon (top) and oxygen (bottom).The dashed horizontal lines indicate the ±10% variations, chosen to correspond to the 1σ uncertainty for these parameters.The interactions were generated using the T2K flux.
of the SRC pair.It requires the energy of the pair to be higher than 14 MeV, i.e. approximately twice the average nucleon removal energy.As a consequence, the SF model implementation in NuWro predicts a larger SRC contribution, amounting to ∼15% of the total CCQE interactions.
The impact of these different implementation choices in NEUT and NuWro on the missing-momentum distribution is displayed in Fig. 9.While the high p m tail in NEUT is exclusively due to SRC, it is not necessarily the case in NuWro because of the additional condition on the energy of the SRC pair.These clear differences motivate the need for a large uncertainty on the SRC contribution to the cross section.This is applied as a normalization parameter of SRC events.We set its prior uncertainty to a value of 100%.Although it should be noted that this uncertainty does not cover the full p m shape differences in Fig. 9. Improved SRC uncertainties would require a more theory-driven SRC implementation which contributes significantly across the full range of p m (to allow phase space coverage for event reweighting).It should also be noted that the SRC alteration parameter allows modification of the total CCQE cross section.There is no attempt in this analysis to offset alterations to the SRCs by modifications to the MF-region shell normalizations.PWIA the model is built on.The NEUT implementation of PB, as described in Eq. 5, provides a simple freedom to vary its treatment, whilst additional uncertainties on the alteration to the cross-section at low energy transfers due to FSI effects can also be considered.

Pauli blocking
As discussed in Sec.II, the SF model in NEUT features a simple description of Pauli blocking (PB) inspired by the RFG model of the nuclear ground state [9].With this approach, events with an outgoing primary nucleon with momentum below the Fermi surface momentum, p F , are removed and do not contribute to the total cross section.In NEUT, p F is set to 209 MeV/c for both carbon and oxygen.PB both reduces the cross section predicted by the SF model and causes significant shape changes for events with low momentum outgoing nucleons, often corresponding to low momentum transfer.Analyses of p F from theory and electron scattering data span a range of ∼20 MeV/c [27,[35][36][37], although we opt for a more conservative size of uncertainty, given the simplistic approach of the PB discussed earlier.Consequently, a parameter varying the threshold p F is prescribed separately for each nuclear target with a conservative prior uncertainty of ±30 MeV/c.Note that NEUT simulates Pauli blocking effects both as part of the primary neutrino interaction and then again as part of its intranuclear FSI cascade simulation.The parameter introduced here only affects the former.Uncertainties related to variations of the FSI cascade are discussed in Sec.III C 2. Fig. 10 shows the impact of varying the PB threshold parameter on the lepton kinematics.Since PB only impacts events with outgoing nucleons with low momentum, its impact is most prominent at low energy trans-fers.Therefore, PB is most noticeable for forward-going leptons.The impact of varying PB is not significant in semi-inclusive measurements because low momentum nucleons are seldom detected due to the high thresholds in current detectors (p p ∼ 450 MeV/c).

Optical potential
The addition of an intranuclear cascade to simulate FSI distorts the outgoing nucleon momentum distribution and accounts for additional hadron ejection, but does not change the inclusive CCQE cross section.A full treatment of the distortion of the outgoing nucleon wave-function would affect the cross section as a function of both lepton and nucleon kinematics [28].To account for this missing alteration to the inclusive cross section due to FSI-which is physics beyond the PWIA-Ref.[27] calculates a correction based on an optical potential (OP) that is derived as a correction to the SF model's prediction of the lepton kinematics.
Although such effects are not implemented in NEUT, NuWro v19.02.01 has an option to include this correction.NuWro was thus used to define an OP parameter that varies the cross section in energy and momentum transfer (q 0 , |⃗ q|).The parameter is the amount of the OP effect on the nominal NEUT MC, from 0% (no correction) to 100% (full correction), using a linear interpolation between the two.The calculation is only available for carbon, although it is not expected to be dramatically different for oxygen, so the parameter is applied independently for each target.To avoid fit convergence issues at the boundaries of this parameter, the prior central value is set to 50% with a ±50% uncertainty.Hence, the prior acts to stabilize the fit whilst being conservative enough that any constraint in the fit is dominated by those im- Event rate [a.u.] FIG. 10.The distributions of the angle between the outgoing muon and neutrino (top) and the muon momentum in the cos θµ > 0.9 region (bottom), showing the impact of varying the PB threshold from −1.5σ (cyan) to +1.5σ (pink), compared to the nominal (black).The interactions were generated on a carbon target using the T2K flux.
plied by cross-section measurement.
As discussed in Sec.II, effects beyond the PWIA appear at low energy and momentum transfer, where the bound nucleons can no longer be treated as independent entities.Similarly to PB, Fig. 11 illustrates the impact of applying the OP correction on the lepton kinematics.The largest impact is again for forward-going leptons, since the transferred momentum and energy in such interactions is the smallest.Whilst the single-dimensional projections show the impact of OP and PB parameters to be similar, they are not degenerate in higher dimensions: FIG.11.The distributions of the angle between the outgoing muon and neutrino (top) and the muon momentum in the cos θµ > 0.9 region (bottom), showing the impact of varying the PB threshold from from 0%, i.e. nominal, (black) to 100% (orange), compared to the nominal (black).
for example, their impact on muon momentum for fixed ranges of muon angle are different.
Since the impact of this correction is calculated only for the outgoing lepton kinematics, the method of applying it to the NEUT distribution of (q 0 , |⃗ q|) from NuWro only alters kinematics of the outgoing nucleon via the leptonnucleon correlations in the implemented SF model.A full treatment of an OP correction would alter these correlations, and so it is expected that the SF model with the OP correction may not have strong predictive power for outgoing nucleon kinematics.An additional uncertainty on the nucleon kinematics is implemented via un-certainties in the FSI cascade is applied, which allows some freedom to mitigate this issue and is discussed in Sec.III C 2.

C. Additional CC0π uncertainties
To benchmark and validate the uncertainty parametrization of the SF model described in the previous sections, cross-section measurements from both the T2K and MINERvA collaborations are fit, as is discussed in Sec.IV.The topology of interest is the charged-current interactions without pions in the final state, known as CC0π, since it provides a sample enriched with CCQE events, thus relevant for studying the SF model.This topology also contains significant fractions of 2p2h and single-pion production events with a subsequent absorption of the outgoing pion, which thereby necessitate additional systematic uncertainties to compare to the measurements.The non-CCQE contributions relevant to the T2K and MINERvA measurements are different due to their different neutrino energy range.T2K has a narrow-band neutrino energy spectrum with a peak neutrino energy of E ν ∼ 0.6 GeV, and MINERvA has a wide-band beam with average neutrino energy ⟨E ν ⟩ ∼ 3.5 GeV.

CCQE axial form factor
Since the treatment of the axial form factor can significantly impact the CCQE cross section-especially at high four-momentum transfers-we also consider the systematic uncertainty related to the dipole form of the axial form factor, typically used in neutrino event generators.The value of the nucleon axial mass M QE A in the dipole form factor is set to 1.03 ± 0.06 GeV/c 2 , estimated from analysis of bubble chamber data [10].Additionally, the dipole approximation of the form factor may be insufficient to describe the evolution of the cross section, particularly at high Q 2 [38].Therefore, three "high Q 2 parameters", which vary the normalization of CCQE interactions in the ranges and [1.00 − ∞] in GeV 2 , are applied to allow freedom beyond the dipole form factors.This approach is derived from the Q 2 > 0.25 GeV 2 parameters in Ref. [10].As also detailed in Ref. [10], the corresponding uncertainties are derived from comparisons of the shape in the high Q 2 region of the dipole and other models, including the z-expansion model [39,40].

Nucleon final-state interactions
As discussed, FSI play a crucial role in altering the outgoing nucleon kinematics and distorting interaction topologies.This is especially relevant when the crosssection measurements are performed in variables that use both the lepton and the hadron information, such as δp T .To account for nucleon FSI uncertainty in this work, CC0π events are divided into two classes: • "With FSI": when the outgoing nucleon kinematics are modified due to FSI, • "Without FSI": in the opposite case, when the nucleon exists the nucleus without being impacted by the intranuclear cascade.
A normalization parameter for each class is applied with a broad 30% prior uncertainty, motivated by an analysis of nuclear transparency measurements [41].The two parameters are very strongly anti-correlated to ensure that the total cross section remains almost constant for each interaction mode.If the amount of "Without FSI" events in a certain interaction mode is reduced, the amount of "With FSI" events in that same interaction mode is therefore increased, yielding a shape-only effect.Effectively these two parameters act almost as a single degree of freedom and are only split for technical implementation reasons, which has a similar impact as changing the mean free path of the nucleon within the nucleus.It should be noted that this approach does inadvertently allow some small changes in shape to the cross section as a function of lepton kinematics, which should remain unchanged by FSI cascades.This is an extremely simple parametrization of FSI uncertainties which is unlikely to be sufficient for a detailed analysis of outgoing hadron kinematics, but can be suitable for analyses that are only sensitive to the coarse hadron kinematic information considered in this work.Improvement of FSI cascade uncertainty parametrizations is a crucially important topic for modeling neutrino interactions in active development (see e.g.Refs [28,42,43]), but is beyond the scope of this work.

Resonant pion production
A non-negligible fraction of CC resonant production (CCRES) contributes to the CC0π topology, which occurs when the outgoing pion is absorbed inside the nucleus through FSI.This is particularly important for multi-GeV neutrinos like in the MINERvA flux, leading the CCRES contribution to be larger than in T2K.To account for this, a normalization parameter that varies the amount of CCRES events with an absorbed pion in this specific topology is considered.Since the impact of FSI is nuclear target dependent, this parameter is split for interactions on carbon and oxygen targets.
In addition, we further consider three parameters that modify the Rein-Sehgal model with lepton mass effects and updated form factors [44][45][46][47][48] implemented in NEUT.The parameters of this model are: the axial mass, M RES A , the value of the axial form factor when Q 2 = 0, C A 5 , and the normalization of the non-resonant 1/2-isospin background, I 1/2 .Their prior values and uncertainties are fixed from recently-updated fits to bubble chamber data on hydrogen and deuterium from ANL [49] and BNL [50], which has been used in Ref. [51,52].

2p2h interactions
As already mentioned, 2p2h interactions also populate the CC0π topology.NEUT describes 2p2h interactions using the Nieves et al. model [53], whose cross section features two distinct peaks in the energy and momentum transfer that broadly correspond to ∆ and non-∆ excitations.
An uncertainty on the amount of 2p2h events is added as a simple normalization parameter which is able to adjust the number of 2p2h interactions for each target.For this parameter, we opt for a 30% prior uncertainty to cover differences between the NEUT model and the alternative 2p2h calculation from the SuSAv2 theory group [22,54].Additionally, we prescribe a 2p2h shape uncertainty which varies the relative strength between the ∆ and non-∆ contributions, as used and described in Refs.[10,55,56].

IV. FITS TO CC0π CROSS-SECTION MEASUREMENTS
Fits to a variety of available neutrino cross-section measurements were performed with the goal of benchmarking the new proposed parametrization; understanding the extent to which these systematic uncertainties can be constrained by data; and assessing the compatibility of these constraints with expectations from electron scattering data.The topology of interest is CC0π, since it is the most sensitive to the CCQE contribution and thus to this new set of parameters.Since most of these parameters can alter both the lepton and nucleon kinematics for CCQE interactions consistently-with the exception of OP and FSI parameters-we fit to crosssection measurements reported in muon kinematics and the transverse kinematic imbalance between the outgoing lepton and nucleons, such as δp T .In addition, since the parameters have been developed both for carbon and oxygen targets, it is particularly interesting to validate them against available data for both elements.For these reasons, the following measurements have been chosen for the fits: the T2K CC0π cross section in muon kinematics (p µ , cos θ µ ) on oxygen and carbon [57], the T2K CC0π cross section in δp T , muon, and proton kinematics on CH [11], and the MINERvA CC0π cross section in δp T on CH [33].The oxygen and carbon measurement from T2K is a simultaneous measurement of the cross section on both targets, and as such provides uncertainties correlated across them.
Each cross-section measurement is not sensitive to all the parametrized uncertainties introduced in Sec.III, and as such not every fit includes all the parameters.Tab.II states when parameters are only included in a sub-set of the fits.The oxygen uncertainties are only fit when considering the T2K measurement that includes an oxygen target.The parameters that affect the low energy transfer region (OP and PB) are only included in cross-section measurements that are sensitive to low energy transfer interactions (i.e.not those that require the observation of a proton in the final state -such that its momentum must be above the 450 MeV/c or 500 MeV/c experimental tracking thresholds (the exact value depends on which measurements are being considered) -as this implies the need for a sufficiently high energy transfer that the implemented corrections to PWIA are not relevant).The parameters that affect the FSI intranuclear cascade are only included in fits sensitive to outgoing nucleon kinematics (i.e.those sensitive to the semi-inclusive CCQE cross section).
Fits are performed using the NUISANCE framework [58].Each of the cross-section measurements used in this work were already implemented in NUISANCE, including the bin-to-bin covariance matrices as published by the experiments.The minimization package used in the fits is MINUIT [59].

A. Chi-squared test-statistic
To perform the desired fits, we developed and implemented in NUISANCE a specific χ 2 test-statistic, with the explicit goal of mitigating the impact of Peelle's Pertinent Puzzle (PPP) [16,17], that is well known to cause a preference for artificially low normalization fit results when a good fit to highly correlated data cannot be found.
Usually, a fit can be made by comparing varied simulated predictions to a cross-section measurement using the measurement histogram B = {B 1 , ..., B n } and covariance matrix M = Cov [{B i }] as provided by the experiments and minimizing a standard test-statistic like the chi-squared: where B MC = (B MC 1 , ..., B MC n ) corresponds to the bins in the histogram of simulated predictions.When the bins in the measurement are highly correlated, fits can converge to simulated distributions with an unphysical low overall normalization when a good fit is not possible, which is known as PPP.Whilst PPP is ultimately just a consequence of correlations in the covariance matrix, the preference for low normalizations partially stems from the χ 2 in Eq. 9 assuming that the absolute uncertainty on each bin of the measurement is independent of its normalization.This implies that the relative uncertainty is larger when fitting to models that predict lower normalizations.When a fitted parametrization cannot match the shape of a measured cross section, the χ 2 will therefore tend to be less poor at low normalizations.If the Gaussian assumption implicit in the standard covariance matrix is correct, then this preference for models with low normalizations is genuine and there is no puzzle (beyond why the fitted parametrization is insufficient).However, it is expected that some types of experimental uncertainty that control the overall normalization of the measured cross section (for example flux or background subtraction systematic uncertainties) should not follow a Gaussian, but rather a log-Gaussian distribution [16].In the case of non-Gaussian uncertainties the χ 2 in Eq. 9 may not represent a good test statistic in a fit.
Methods to work around PPP have been employed, such as using a shape-only chi-square [12,31], or neglect-ing the bin-to-bin correlations [60].However, neither of these methods are satisfactory since the former ignores the valuable information on the total cross section from the measurement, while the latter is almost meaningless when the measured cross section has significant correlations between bins (as most do).We propose an alternative way to mitigate this effect which consists of separating the normalization and shape contributions of the covariance matrix, that corresponds to isolating the relative uncertainty from the uncertainty on the overall cross-section normalization, in such a way to make the absolute uncertainty larger for models predicting lower normalizations.This method results in a relative uncer-tainty that remains constant as a function of the normalization, as motivated by the arguments of Ref. [16].Such a treatment is generally well motivated for data dominated by correlations due to multiplicative uncertainties, such as the overall normalization uncertainties from flux uncertainties that often dominate cross-section measurements.
Concretely, this can be obtained by applying a transformation to both the data and simulated histograms as well as to the covariance matrix, in order to separate them into a "shape" and a "norm" part.The new histograms C = {C 1 , ..., C n } are defined as: where α is a scale parameter.Note that if B is a differential rather than absolute number of events, scaling by the bin width is required.The function f : B → C is bijective, meaning that no information is lost by moving from B to C.
In this new basis, a covariance matrix, N = Cov [{C i }], would ideally be built using the same "toys" or "universes" experiments typically use to build their standard covariance.Unfortunately such information is not usually provided in experimental data releases.Instead, the covariance matrix is approximated via a non-linear transformation of the original covariance matrix M , using the following formula: where J(f ) is the Jacobian of the non-linear transformation f .The new covariance matrix is expressed as follows: The matrix N has the same dimension and the same positive-definiteness properties as M since the mapping B → C is a bijection.N is composed of two diagonal blocks: the N S block which corresponds to the shape-only covariance, and the kl M k,l element which corresponds to the data norm variance.The off-diagonal blocks represent the correlations between the norm and the shape components.Finally, by transforming the MC histogram using this same function f , the "norm-shape" (NS) chi-square can be computed in this basis as: This new computation of the covariance matrix and the chi-square was implemented in NUISANCE and the correctness of the implementation was validated by comparisons with covariance matrices computed from toys.The use of the χ 2 NS allows to mitigate the PPP problem observed when using the standard χ 2 of Eq. 9 in the fit.This method was first discussed in Ref. [61,62] and also used, via our implementation, in Refs.[60,63] 3 .
3 On a historical note: the MiniBooNE collaboration had de-It should be noted that whilst this method mitigates the impact of PPP, it also makes assumptions that: • the real uncertainty in the data did not follow the multivariate Gaussian covariance reported by the experiments, but rather follows a distribution from which the relative uncertainty is constant as a function of the measurement's normalization, • the transformation of the experimentally reported covariance provides a better description of the real distribution of the measurement uncertainties.
Given the limited information provided in experimental data releases, it is not possible to test these assumptions 4 .As such, the modified test-statistic is used for fits but the usual χ 2 from Eq. 9 is also reported.
veloped (but not published) a similar approach in the context of decomposing covariance matrices into total, shape-only and "mixed" terms.This is noted on page 217 of Ref. [64]. 4 If experiments were to provide a distribution of cross-section results mapped by their uncertainties (i.e. the "universes" used to calculate the covariance matrix that is usually supplied) it may be possible to tailor test statistics used in fits to best describe the plausible variations of the measured cross section.

B. Fit results
In this section, the results of the fits, using the χ 2 NS defined in the previous section, to T2K and MINERvA data are presented and discussed.
1. Fit to T2K CC0π cross section data on oxygen and carbon in muon kinematics Fig. 12 shows the prefit and postfit cross-section predictions for the T2K CC0π measurement on carbon and oxygen targets from Ref. [57], whilst Fig. 13 and Fig. 14 show the prefit and postfit parameter values and the associated correlation matrix.The prefit, i.e. the nominal spectra predicted by NEUT, and the postfit distributions are displayed and compared with the data and the value of the corresponding χ 2 NS as well as the usual χ 2 are reported in Tab.III.
It can immediately be noted that the data/model agreement is dramatically improved after the fit adjustment of the systematic parameters, which is also reflected in the important decrease in both χ 2 NS and χ 2 .Importantly, both are also lower than the number of bins, implying a quantitatively good agreement between the postfit model and the measured cross section.It should be noticed that, as detailed in Ref. [57], the prefit data/model disagreement comes mainly from the bins that correspond to forward lepton kinematics (i.e. the most forward cosθ µ slice).This is the region that corresponds to a low energy transfer which is known to be more complicated to model due to the effects beyond the PWIA.For instance RFG and LFG models are often corrected using the random phase approximation to achieve better agreement in this region of kinematic phase space [29].It is thus expected to see that the postfit agreement in this region is driven by the parametrization of physics beyond PWIA discussed in Sec.III B: PB and OP.Indeed, Fig. 13 shows the preference for an application of the OP correction and more PB than is in the nominal NEUT.Interestingly, the strength of these corrections to PWIA appear to be stronger on carbon than on oxygen.
Fig. 13 shows that the measurement offers fair sensitivity to constraining the shell normalization parameters, which assist in compensating the effects of the PB and OP parameters.The overall CCQE normalization and its shape in Q 2 is also adjusted by a significant variation of the M QE A and high Q 2 parameters.Almost all parameters remain reasonable given their prior expectations both as encoded in the prefit uncertainties and as would be expected from electron scattering data.In particular, the shell normalizations are not pulled far from their nominal values.However, it can be noted that PB for carbon is pulled outside of its relatively conservative prior value to give an effective Fermi momentum of ∼250 MeV/c.Interpreted as a Fermi-gas Fermi momentum, this would imply nonphysical nuclear density [65] and so this is may be indicative that PB is acting partially as an effective parameter to account for missing freedom within the model (e.g. to span differences between different approaches to PB or other physics beyond PWIA).
The postfit correlation matrix is shown in Fig. 14.There are anticorrelations between the PB and OP uncertainties for carbon and more prominently for oxygen since both parameters have a similar impact on the cross section.Anticorrelations are also been between the pion absorption parameter and the parameters that control the CCQE normalization, indicating that the data is not able to offer isolated constraints on the QE and non-QE model components (which is to be expected when only measuring outgoing lepton kinematics).
It may be noteworthy to highlight that the parametrization developed in this work treats the oxygen and carbon uncertainties as two independent groups (except for the nucleon-level parameters such as M QE A , M RES A , etc., which are common), but the postfit correlation matrix shown in Fig. 14 exhibits correlations between these two sets of parameters.As previously stated, the measurement corresponds to a joint cross section on the two targets, the data covariance includes correlations between carbon bins and oxygen bins, which is then reflected on the parameters.Such simultaneous measurements are increasingly important for oscillation measurements to better understand the extrapolation of model constraints between targets.

Fit to T2K cross section data in CC0π0p and CC0πN p topologies on hydrocarbon
In this section, we fit T2K CC0π cross-section measurements on hydrocarbon as a function of transverse momentum imbalance δp T and final-state muon kinematics for CC0π topologies with (CC0πN p, N ≥ 1) and without (CC0π0p) a measured proton in the final state respectively [11].Note that 0p and N p refer to protons above and below tracking threshold (∼500 MeV/c).Three distinct fits are performed: One of the interests of performing this simultaneous fit is to evaluate the ability of the model to describe neutrino interactions in different neutrino energy ranges and in different regions of lepton kinematics.Indeed, the CC0π0p and the CC0πN p topologies correspond to distinct regions of momentum transfer (as a proton requires ∼450 MeV momentum to be observed) and therefore to different regions of E ν , as illustrated in Fig. 15.It also tests the parametrization's ability to describe how     III.Summary of the prefit and postfit norm-shape chi-square χ 2 NS used in the minimization, and the usual chi-square χ 2 for reference in the different fits presented in Sec.IV B along with the corresponding number of bins.
lepton kinematics changes as a function of the outgoing proton kinematics.It should be noted that the two measurements considered in (c) are treated as uncorrelated (as correlations between the two measurements are not available) despite not being completely independent.The kinematic overlap between them is expected to be very small, given that CC0π0p topology allows protons up to 500 MeV/c, whilst the δp T measurement considers only protons above 450 MeV/c although a more complete analysis would provide correlations from systematic uncertainties.
In our analysis the highest outgoing muon momentum bin (which extends up to 30 GeV/c) in each angular slice of the lepton kinematics measurement in CC0π0p is removed.These constitute a negligible fraction of the total measured cross section and are in a momentum range where T2K cannot guarantee reliable reconstruction 5 .
Figs. 16 and 17 show the prefit and postfit distributions from fits (a) and (b) respectively, whereas Figs.18 and 19 display the prefit vs. postfit parameters and correlations respectively from the three fits.As previously discussed and shown in Tab.II, it should first be noted that the three fits do not use the same set of parameters.
The binning of δp T , as shown in Fig. 17, is rather coarse, and so even the prefit chi-squares reported in Tab.III for the CC0πN p fit are relatively low.Nevertheless, the uncertainty parametrization presented in this work yields a significantly improved postfit agreement, as demonstrated by the chi-squares.This is achieved largely thanks to the SF model shell parameters as indicated in Fig. 18 (black).This fit also shows sensitivity to SRC, nucleon FSI and 2p2h shape parameters, as seen by the reduction of their postfit uncertainties.Indeed, these are effects that are probed by δp T , particularly in the tail of its distribution [32].The correlations between them that appear in the postfit correlation matrix (top left of Fig. 19) indicate that these effects cannot be disentangled solely with δp T .The T2K measurement considered here does not attempt to separate the 2p2h and FSI contributions to the tail of the distribution.There is scope to achieve a better separation by performing a measurement of δp T as a function of δα T (another transverse variable [32]), or even as a function of outgoing lepton kinematics.Such a measurement has recently been performed by the MicroBooNE collaboration [34], and there are prospects for similar high-precision measurements using the upgraded T2K ND280 detector [66].On the other hand, even with the coarse binning, the CCQE dominance in the bulk of δp T (which, for CCQE events is just the transverse projection of the missing momentum) allows a notable constraint on the shell normalization and shape parameters.
Contrary to the CC0πN p fit and the oxygen + carbon analysis from Sec. IV B 1, fits including the CC0π0p measurement result in relatively poor chi-squares.This indicates that the model presented in this work is sufficient for analyzing results integrated over outgoing proton kinematics or binned coarsely in δp T , but is insufficient to fully describe the CC0π0p measurement that includes only protons below 500 MeV/c in its signal definition.This is likely due to insufficient freedom to alter the ratio of events with and without a proton above 500 MeV/c as a function of the outgoing lepton kinematics.It is likely that the lack of detailed FSI uncertainties are at least part of the cause of this.
Similarly to the fit of T2K data on oxygen and carbon discussed in Sec.IV B 1, the improved agreement of the model with the data in the CC0π0p measurement is largely driven by the increase in the PB parameter, as shown in Fig. 18 (blue), which affects the forward angular region where the discrepancies are the largest.However, in this case the PB parameter remains within the prior uncertainty range.It can also be seen that OP actually moves to apply a weaker suppression of the low energy transfer region.Overall this seems to suggest that the CC0π0p measurement sees less of a requirement of a forward lepton angle reduction of the cross section.This might suggest that the apparent required CCQE suppression from the fit integrated over nucleon kinematics may have been acting in lieu of a need to reduce the non-QE contributions (which should be reduced in the CC0π0pN p result) in the forward angle region.
Fig. 18 additionally shows that the Q 2 parameter that affects the region 0.25 ≤ Q 2 < 0.50 GeV 2 converges to a value significantly outside of its prior uncertainty when considering the CC0π0p measurement.Note that a smaller raising of the parameter was seen the carbonoxygen measurement considered in Sec.IV B 1. Overall this might indicate a lack of freedom to vary the cross section at high momentum transfer.
The postfit values of the parameters from the simultaneous fit of CC0π0p and CC0πN p measurements are also displayed in Fig. 18 (purple).Most of the parameters converge to similar values as in the CC0π0p-only fit since the corresponding data statistically dominates the chi-square and drives the fit.In particular, the FSI parameters are less pulled by the fit in comparison with the CC0π0p-only one thanks to the constraint from the δp T distribution tail.The bottom panel of Fig. 19 shows strong (anti)correlations between parameters related to PB, OP, FSI and multinucleon effects.This explains the slightly different postfit values for some of these parameters between the CC0π0p-only and the simultaneous fit.With a more statistically significant δp T measurement, we could expect a better separation between the struck nucleon-related uncertainty (FSI, SRC, 2p2h) probed by δp T and the low energy transfer effects probed mostly by the forward outgoing muon kinematics.

Fit to MINERvA cross section data in CC0πN p topology on hydrocarbon
Results from the fit to the MINERvA measurement of the CC0πN p (considering N protons above 500 MeV/c) cross section as a function of δp T are shown in Fig. 20.As suggested by the chi-square values quoted in Tab.III, the data-MC prefit agreement is quite poor and is slightly improved in the postfit.Nevertheless, Fig. 21 exhibits a clear sensitivity to most of the considered parameters, including the missing-momentum shape uncertainties.This is partially due to the significantly finer binning in the MINERvA data in comparison with the T2K measurement.This allows a more precise probe of the nuclear effects that impact the δp T distribution.In fact, as discussed in Sec.II, the bulk is sensitive to Fermi motion which is mainly affected by the shell normalization and shape parameters (see Figs. 4 and 8).On the other hand, its tail can be altered by SRC, 2p2h, CCRES and FSI uncertainties.
However, it is clear from the relatively high value of the postfit χ 2 NS that the present parametrization of the SF FIG.18. Prefit and postfit values and constraints on the uncertainties from the fit to T2K CC0π0p measurement of lepton kinematics and CC0πN p measurement of δpT on carbon.The displayed central value for each parameter corresponds to the difference with respect to its prior value divided by the prior uncertainty as reported in Tab.II.systematic uncertainties does not provide enough freedom in the model to entirely cover discrepancies with the measurement.This can be attributed to the fact that, due to the higher energy of the MINERvA flux, there is a significant contribution from the other interaction channels (like CCRES) through pion absorption that would need a more sophisticated parametrization.For instance, the predicted fraction of CCRES events by NEUT corresponds to almost ∼ 20% of the CC0π topology for MINERvA, which is below 10% in the case of T2K.This larger CCRES component is also responsible for the tighter constraints on the Rein-Sehgal parameters in comparison with the previous fits.The relatively poor agreement can also mean that the current parametrization of the CCQE model may need further improvements, especially for FSI effects.
Fig. 22 reports the postfit correlation matrix for this fit.In comparison with the T2K δp T fit, we can notice that the anticorrelation between the p-and the s-shell normalization parameters is less prominent thanks to the finer binning which alleviates their degeneracy with a more precise probe of the shape of the δp T distribution.Besides, the SRC, 2p2h, FSI and CCRES uncertainties are correlated as expected since they affect the same highδp T region.

V. IMPLICATIONS FOR OSCILLATION ANALYSES
In order to qualitatively evaluate the impact of this new parametrization of uncertainties for the SF model in the context of future neutrino oscillation measurements, we can consider the prefit and postfit uncertainties projected onto observables oscillation analyses are particularly sensitive to: namely the neutrino energy dependence of the cross section and the E QE ν bias.Fig. 23 shows the prefit and postfit spectra and con-straints for the distribution of true neutrino energy and the E QE ν bias expected in the T2K flux.These are split into showing the total constraint on a differential cross section (left) and on only its shape (right).This split is informative as it allows a separation of the overall constraint placed on the total normalization of the cross section, which should be relatively independent of the uncertainty parametrization used, from the constraint placed on the shape of the distributions.The postfit (prefit) distributions are obtained using an ensemble of 500 distributions sampled from the posfit (prefit) values and covariance from the fit shown in Sec.IV B 1.
The postfit constraints on the cross section as a function of E ν are significantly improved in comparison with the prefit, as shown in the left plots of Fig. 23.This is more visible in the bottom-left panel for the bias of E QE ν particularly around the true E ν (i.e.0 on the plot), which is the region that is most affected by CCQE events and the SF model.The negative tail, which is more affected by multinucleon effects and CCRES interactions, is only slightly impacted since the measurement used in the fit has only a small component of these interactions.
It is clear from the bottom-right plot of Fig. 23 that the uncertainty model offers significant freedom in the shape of the neutrino energy bias and that this is well constrained from the fit to the T2K cross-section measurement.On the other hand, the top-right plot shows that the freedom in the shape of the neutrino energy dependence of the cross section is more limited and is not so strongly constrained by the fit.
To quantify the impact of the reduced uncertainties on the allowed space in the differential cross section as a function of the neutrino energy (top left panel of Fig. 23), we can evaluate the determinant of the prefit and the postfit bin-to-bin covariance matrices of n b dimensions, where n b = 10 is the number of E ν bins.In fact, the square root of the determinant of the covariance is proportional to the volume of the n b -dimensional ellipsoid  that covers the N -sigma variations of the Gaussian errors.We find that this volume is reduced by a factor of 4.68 thanks to the constraints from the fit.

VI. PERSPECTIVES AND CONCLUSIONS
The need to use a more sophisticated parametrization of nuclear model uncertainties in neutrino oscillation experiments is becoming increasingly urgent.In this paper, we introduce a parametrization of systematic uncertainties on the inclusive and semi-inclusive predictions of the Benhar spectral function model driven mostly by the natural degrees of freedom within the model, as well as a few additional freedoms to account for the limitation of SF's PWIA description of neutrino interactions.Fits to existing T2K and MINERvA CC0π cross-section mea-surements show that this parametrization is able to offer a well-motivated means to improve the agreement with respect to the nominal model predictions, especially at the T2K energies.The fits were able to achieve a quantitatively good postfit agreement with T2K measurements of outgoing lepton kinematics on carbon and oxygen targets as well as for T2K measurements of δp T , but the postfit agreement is less satisfactory when fitting T2K measurements of lepton kinematics within a restricted proton kinematic phase space.Overall this suggests that uncertainty parametrization is likely to be broadly sufficient for describing the inclusive CCQE cross section at T2K energies (although care is required with regards to the treatment of the parameters governing the cross section at low energy transfer), but that additional components are likely to be needed when detailed modeling of hadron kinematics is required.The SF model has been adopted by the T2K collaboration for neutrino oscillation analyses since 2020 and the described parametrization has been used for the latest iteration of its oscillation analysis [51,52], which focuses on measurements of primarily outgoing lepton kinematics.In the longer term, future improved measurements at T2K will be possible thanks to the new detectors like the Super-FGD in the upgrade of the T2K near detector [6].Its fine-grained design will allow to precisely measure the kinematics of the hadronic products from neutrino interactions, allowing for instance the reconstruction of protons with momenta down to 300 MeV/c.Ref. [15] shows quantitatively the expected improvements of the constraints on the nuclear effects described in the SF model with a simplified version of the parametrization introduced in this paper.It demonstrates how exploiting nucleon-lepton correlations can considerably constrain the uncertainties discussed in this work thanks to not only measurements of the single-transverse variables, but also an improved estimator of neutrino energy based on the sum of muon energy and nucleon kinetic energy.
Beyond T2K, the proposed model parametrization is likely to form a reasonable starting point for other experiments that use nuclear shell models.This includes experiments that measure exclusive final states, as most of the added parameters offer a means to alter both lepton and hadron kinematics in a consistent way, although semi-inclusive cross-section measurements suggest some extensions to the parametrization, in particular related to freedoms within the FSI model, will likely be required.(bottom) differential cross sections in the T2K beam on a carbon target following the fit to the T2K CC0π joint measurement of lepton kinematics on carbon and oxygen (Sec.IV B 1).The plots on the left show the overall constraint on the differential cross section, whilst the plots on the right indicate the constraint on the shape of the distribution.

FIG. 1 .
FIG.1.Simulations of the two-dimensional distribution of missing energy, Em, and missing momentum, pm, for carbon using the RFG (blue), LFG (green) and SF (red) model simulated with the NEUT neutrino interaction generator.

FIG. 2 .FIG. 3 .
FIG. 2.Comparison between the distribution of the missing momentum (left) and the neutrino energy bias (right) for CCQE interactions for the RFG (blue), LFG (green) and SF (red) model generated by NEUT using the T2K muon neutrino flux.

FIG. 4 .
FIG.4.The distributions for removal energy (top) and missing transverse momentum (bottom) as given by the SF model implemented in NEUT, showing the impact of the shell normalization parameters (N shell ) from −1.5σ (pink) to +1.5σ (cyan), compared to the nominal (black).The interactions were generated on a carbon target using the T2K flux.

FIG. 5 .
FIG.5.Distribution of the missing momentum within each nuclear shell for carbon as given by the SF model implemented in NEUT generated using the T2K flux.

FIG. 7 .
FIG.7.Distributions of the missing momentum from the NEUT SF inputs (black) compared to electron scattering measurements[20] (blue and red), made for different nucleon and lepton kinematics for carbon in the p-shell (left) and the s-shell (right).

FIG. 9 .
FIG. 9. Distribution of the missing momentum in the SF model in NEUT (left) and NuWro (right) for carbon with the MF and SRC contributions.
FIG.12.Prefit (red) and postfit (blue) distributions of pµ in bins of cos θµ from the fit to T2K CC0π joint measurement of lepton kinematics on carbon and oxygen.The usual chi-squares as well as the number of bins are quoted in the legend.The NS chi-square χ 2 NS used in the minimization is reported in Tab.III.

FIG. 14 .
FIG.14.Postfit correlation matrix from the fit to T2K CC0π joint measurement of lepton kinematics on carbon and oxygen.

FIG. 15 .FIG. 16 .FIG. 17 .
FIG. 15.True squared four momentum transfer (upper) and neutrino energy (lower) distributions of interactions falling into the the CC0π0p (blue) and CC0πN p (red) topologies in the T2K beam on a carbon target as predicted by NEUT.

FIG. 19 .
FIG.19.Postfit correlation matrices from the fit to T2K CC0π0p measurement of lepton kinematics (top left) and CC0πN p measurement of δpT (top right), as well as the simultaneous fit (bottom).

FIG. 20 .FIG. 21 .
FIG.20.Prefit (red) and postfit (blue) distributions from the fit to MINERvA CC0π + N p measurement of δpT on carbon.The usual chi-squares as well as the number of bins are quoted in the legend.The NS chi-square χ 2 NS used in the minimization is reported in Tab.III.

FIG. 22 .
FIG. 22. Postfit correlation matrix from the fit to MINERvA CC0πN p measurement of δpT on carbon.

FIG. 23 .
FIG. 23.Prefit (red) and postfit (blue) constraints on the true neutrino energy (top) and the bias of E QE ν

TABLE I .
Target Shell E shell [MeV] σ shell [MeV] N shell uncertainty Nuclear energy levels with their widths for the different shells, for the SF as implemented in NEUT.The last column represents the relative prior uncertainty set on the corresponding shell normalization parameter, which all have a central value of 0.

TABLE II .
Summary of the parameters introduced in Sec.III and their prior uncertainties.
13G.13.Prefit (red) and postfit (blue) values and constraints on the uncertainties from the fit to T2K CC0π joint measurement of lepton kinematics on carbon and oxygen.The displayed central value for each parameter corresponds to the difference with respect to its prior value divided by the prior uncertainty as reported in Tab.II.