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

Parametrized uncertainties in the spectral function model of neutrino charged-current quasielastic interactions for oscillation analyses J. Chakrani ,1,* S. Dolan ,2,†M. Buizza Avanzini ,1,‡A. Ershova ,3L. Koch ,4K. McFarland ,5G. D. Megias ,6 L. Munteanu ,2L. Pickering,7K. Skwarczynski ,8V. Q. Nguyen,1and C. Wret 9 1Laboratoire Leprince-Ringuet, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, Palaiseau, France 2European Organization for Nuclear Research (CERN), Geneva, Switzerland 3IRFU, CEA, Universit´e Paris-Saclay, Gif-sur-Yvette, France 4Institut für Physik, Johannes Gutenberg-Universität, Mainz, Germany 5University of Rochester, Department of Physics and Astronomy, Rochester, New York, USA 6Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla, Spain 7Royal Holloway University of London, Department of Physics, Egham, Surrey, United Kingdom 8National Centre for Nuclear Research, Warsaw, Poland 9Oxford University, Oxford, United Kingdom (Received 5 October 2023; accepted 11 March 2024; published 15 April 2024) 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 propose such a set of parameters, which can be applied to nuclear shell models, and test their application to the Benhar et al. [Nucl. Phys. A579, 493 (1994)] 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. DOI: 10.1103/PhysRevD.109.072006 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ν∼1GeV 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 finalstate 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, nonstatic 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 *[email protected] †[email protected] ‡[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 109, 072006 (2024) 2470-0010=2024=109(7)=072006(25) 072006-1 Published by the American Physical Society 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 charged-current 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–14].It has also been used to study the sensitivity of T2K’s neardetector 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 parametrization is benchmarked by fitting to neutrino-nucleus cross section measurements from the T2K and MINERvA experiments with a view to determining which data is able to be quantitatively well-described and to explore the constraints that may be expected from near-detector data in neutrino oscillation analyses. 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. Va 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 neutrino-nucleus 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 “shortrange correlations”(SRC) between nucleons within the nucleus. Figure 1shows a comparison of the predictions between the three models for the distribution of the nucleon removal energy Ermv, and the initial-state nucleon momentum pN, 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→lþN0þA0, where Aand A0are the initialstate nucleus and the final-state remnant respectively and N0is the outgoing nucleon), these can be defined as1 FIG. 1. Simulations of the two-dimensional distribution of the removal energy Ermv, and the initial-state nucleon momentum pN, for carbon using the RFG (blue), LFG (green), and SF (red) model simulated with the NEUT neutrino interaction generator. 1The removal energy and the initial-state nucleon momentum defined here are the same as the missing energy and the missing momentum introduced in Ref. [9]. It should be noted that this equivalence is only accurate for the special case of single nucleon knock out with no FSI. The missing energy and momentum observables reported by electron scattering experiments (as in Ref. [20]) must usually be corrected in order to obtain removal energy and initial-state momentum. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-2 Ermv ¼Eν−El−TN0−TA0−ΔmN;ð1Þ pN¼j  pNj¼j  pν− pl− pN0j;ð2Þ where Eνand  pν(Eland  pl) are the incoming neutrino (outgoing charged lepton) energy and momentum, respectively, TN0and pN0are the kinetic energy and momentum of the pre-FSI outgoing nucleon respectively, TA0is the kinetic energy of the nuclear remnant, and ΔmNaccounts for the mass difference of the initial- and final-state nucleon. TA0is calculated from pN0assuming a mass of the ground state of the remnant nucleus following the interaction. Figure 1clearly highlights the nuclear shell structure of carbon encoded in the SF model, with a sharp p-shell at Ermv ∼18 MeV and a diffuse s-shell at Ermv ∼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 pNis uniformly distributed in 3-momentum,  p, below the Fermi momentum pF, which yields a cliff feature at pN∼220 MeV=c. The LFG and SF models both predict narrower pNdistributions than the RFG model. A highmomentum 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 Eb, and that there is a single outgoing nucleon,2 EQE ν¼2El ˜ mN−ðm2 lþ˜ m2 N−m2 NÞ 2ð˜ mN−Elþplcos θlÞ;ð3Þ where mNis the mass of the struck nucleon, ˜ mN¼ mN−Eb, and ml,El,pl, and θlare 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 ground-state models. The smearing is dominated by the isotropic Fermi motion of the bound nucleon. The simplistic RFG model is offset due to its shifted Ermv 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, EQE ν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 removal energy, Ermv, and the initial-state momentum, pN, of the nucleons within the nucleus. It is implemented in FIG. 2. Comparison between the distribution of the initial-state nucleon 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. 2In the case of single-track events intended to isolate single pion events, the outgoing nucleon mass can be replaced by the Δð1232Þmass. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-3 both the NEUT and NuWro event generators, where the joint probability distribution of Ermv and  pNof nucleons within the nucleus is written as Pð pN;E rmvÞ¼PMFð pN;E rmvÞþPcorrð pN;E rmvÞ;ð4Þ where PMFð pN;E rmvÞis the mean field (MF) part and Pcorrð pN;E rmvÞ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 pN(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  pNto 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 pNto 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 Ermv,pNspace. The two-dimensional (pN,Ermv) 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 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 Pð pN;E rmvÞ→Pð pN;E rmvÞθðj pNþ qj−pFÞ;ð5Þ where θis the Heaviside step function,  qis the transferred momentum, and pFis 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. 0 50 100 150 200 250 300 350 400 450 [MeV/c] N p 0 50 100 [MeV] rmv E 0 50 100 150 200 250 300 350 400 450 [MeV/c] N p 0 50 100 [MeV] rmv E FIG. 3. The prediction of the two-dimensional distribution of removal energy and initial-state nucleon momentum for carbon (left) and oxygen (right) for NEUT. The brightness of the color represents the probability of finding an initial-state nucleon with a particular removal energy and momentum state. The white lines indicate the cuts used to separate the MF region (low Ermv,pN) from the SRC region (high Ermv,pN) in the NuWro (dashed) and NEUT (solid) neutrino interaction generators. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-4 PB is a first-order effect beyond the impulse approximation part of PWIA. Whilst PWIA is justifiable for highenergy transfer, the interactions with low-energy transfer (q0≲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 reinteraction 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 favor 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 single-transverse variables [32] allow a probe of nuclear effects, particularly with the measurement of the transverse momentum imbalance δpTdefined from the kinematics of the outgoing charged lepton and nucleon as δ pT¼ pT lþ pT N0;ð6Þ where  pT land  pT N0are 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 δpTtherefore 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 δpTshow 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 (e.g., background and final-state interaction corrections made when inferring spectral functions from electron scattering measurements [35]). The exact priors used are stated and justified in the following subsections. Overall, Table 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. 1. Shell occupancy As shown in Fig. 3, the removal 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 differ from those predicted by a simple independent-particle shell model [25,26,35–37]. The observed difference is often explained as an exploration of the constraints inherent in PWIA and in particular the PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-5 significance of short- and long-range nucleon-nucleon correlations. This variation is commonly assessed using spectroscopic factors which gauge the actual occupancy of each shell. Spectroscopic factors obtained from electron scattering data generally fall within the range of 0.5–0.7, contingent on the specific orbital and nucleus under consideration [38]. Whilst as effective quenching is naturally included within the SF model (being derived from electron scattering data), its determination relies on modeldependent corrections to measurements, in particular to account for final-state interactions [35], which come with significant uncertainty. 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: fshellðErmvÞ¼1þNshell × exp −ðErmv −EshellÞ2 2σ2 shell ð7Þ ¼1þNshell ×gshellðErmvÞ;ð8Þ where Nshell is the normalization parameter of a given shell, and Eshell and σshell correspond to the center and the width of the Gaussian function, gshellðErmvÞ, 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 Eshell and σshell are derived from an analysis of the removal energy distributions in NEUT, and are shown for each shell in Table I. The impact of varying Nshell for the carbon p- and s-shells is shown in Fig. 4. As can be seen in the Ermv distribution from Fig. 4, this approach results in a thinner or thicker shell for higher or lower Nshell 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 Nshell). Note also that there is no limit to the range of Ermv 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 Nshell to give the event a significant weighting. One important effect of the shell occupancy parameters is how they alter the total pNdistribution, since the initialstate 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 pN, and consequently the distribution of δpT, 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. FIG. 4. The distributions for removal energy (top) and transverse momentum imbalance (bottom) as given by the SF model implemented in NEUT, showing the impact of the shell normalization parameters (Nshell) 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. TABLE I. 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. Target Shell Eshell (MeV) σshell (MeV) Nshell uncertainty Carbon p18 15 0.2 s36 25 0.4 Oxygen p1=212 8 0.45 p3=219 8 0.25 s42 25 0.75 J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-6 As discussed, we opt to keep the prior uncertainties on the Nshell parameters conservative. The chosen values are reported in Table I. These values cover differences in the shell occupancy of the SF model [derived from ðe; e0pÞ 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%). Figure 6shows the impact of varying each shell normalization parameter Nshell on the total CCQE cross section, and Fig. 4demonstrates potential sensitivity to variations through measurements of δpT 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. 2. Initial-state nucleon momentum shape Electron scattering data from Ref. [20] provides measurements of the initial-state nucleon momentum distributions pNfor carbon. By comparing the predicted NEUT distributions of pNin 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 initial-state nucleon momentum which correspond to ðe; e0pÞkinemat- ics with the most different four-momentum transfer, Q2, compared to the NEUT prediction. This builds an initialstate nucleon momentum shape uncertainty for each shell that changes the shape of the pNdistribution. Each shell parameter is defined between -1 and 1, corresponding to the two extreme pNdistributions, and a linear extrapolation is implemented beyond the range ½−1;1. These pNshape parameters are expected to mainly affect experimental observables sensitive to the initial-state 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 δpTis affected by this uncertainty since it corresponds, in the absence of other nuclear effects, to the transverse projection of the Fermi motion. 3. 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 FIG. 5. Distribution of the initial-state nucleon momentum within each nuclear shell for carbon as given by the SF model implemented in NEUT generated using the T2K flux. FIG. 6. The impact of the shell normalization parameters (Nshell) 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. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-7 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 Ermv and pNto 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 Ermv >100 MeV or pN>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 ðpN;E rmvÞ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 knockout 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 initial-state nucleon momentum FIG. 8. The distributions for the angle between the outgoing lepton and neutrino (left) and the transverse momentum imbalance (right) as given by the SF model implemented in NEUT, showing the impact of variations of the p-shell shape uncertainty from −2σ(pink) to 2σ (cyan), compared to the nominal (black). The interactions were generated on a carbon target using the T2K flux. FIG. 7. Distributions of the initial-state nucleon 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). J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-8 distribution is displayed in Fig. 9. While the high-pNtail 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 pNshape differences in Fig. 9. Improved SRC uncertainties would require a more theory-driven SRC implementation which contributes significantly across the full range of pN(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. B. Uncertainties for physics beyond PWIA In addition to parameters to alter the natural degrees of freedom within the SF model, additional freedoms are required to account for plausible variations from the 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. 1. Pauli blocking As discussed in Sec. II, the SF model in NEUT features a simple description of 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, pF, are removed and do not contribute to the total cross section. In NEUT, pFis 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 lowmomentum outgoing nucleons, often corresponding to low-momentum transfer. Analyses of pFfrom theory and electron scattering data span a range of ∼20 MeV=c [27,39–41], 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 pFis 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. Figure 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 transfers. Therefore, PB is most noticeable for forward-going leptons. The impact of varying PB is not significant in semiinclusive measurements because low-momentum nucleons are seldom detected due to the high thresholds in current detectors (pp∼450 MeV=c). 2. 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 FIG. 9. Distribution of the initial-state nucleon momentum in the SF model in NEUT (left) and NuWro (right) for carbon with the MF and SRC contributions. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-9 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 Table III. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-16 FIG. 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 Table II. FIG. 14. Postfit correlation matrix from the fit to T2K CC0πjoint measurement of lepton kinematics on carbon and oxygen. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-17 (iii) simultaneously fitting lepton kinematics in CC0π0p and δpTin CC0πNp. 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π0pand the CC0πNp 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 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π0ptopology allows protons up to 500 MeV=c, whilst the δpTmeasurement 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π0pis 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 Figures 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 Table II, it should first be noted that the three fits do not use the same set of parameters. The binning of δpT, as shown in Fig. 17, is rather coarse, and so even the prefit chi-squares reported in Table III for the CC0πNp 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 δpT, particularly in the tail of its distribution [32]. The correlations between them that appear in the postfit TABLE III. Summary of the prefit and postfit norm-shape chi-square χ2 NS used in the minimization, and the usual chi-square χ2for reference in the different fits presented in Sec. IV B along with the corresponding number of bins. Measurement Prefit χ2 NS Postfit χ2 NS Prefit χ2Postfit χ2Number of bins T2K oxygen þcarbon (Sec. IV B 1) 110.88 35.81 98.79 30.30 58 T2K CC0πNpδpTonly (Sec. IV B 2) 12.59 7.37 15.72 8.48 8 T2K CC0π0pðpμ;cos θμÞonly (Sec. IV B 2) 144.35 75.13 107.57 62.55 50 T2K CC0π0pþCC0πNp (Sec. IV B 2)144.35 þ14.56 86.80 þ10.01 107.57 þ16.76 64.19 þ11.83 50 þ8 MINERvA δpT(Sec. IV B 3) 109.10 79.51 114.32 76.14 24 FIG. 15. True squared four momentum transfer (upper) and neutrino energy (lower) distributions of interactions falling into the CC0π0p(blue) and CC0πNp (red) topologies in the T2K beam on a carbon target as predicted by NEUT. 5Whilst these bins are in the data release, they were not shown within Ref. [11]. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-18 correlation matrix (top left of Fig. 19) indicate that these effects cannot be disentangled solely with δpT. 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 δpTas 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 [70]. On the other hand, even with the coarse binning, the CCQE dominance in the bulk of δpT(which, for CCQE events is just the transverse projection of the initial-state nucleon momentum) allows a notable constraint on the shell normalization and shape parameters. Contrary to the CC0πNp 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 δpT, but is insufficient to fully describe the CC0π0pmeasurement 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 FIG. 16. Prefit (red) and postfit (blue) distributions of pμin bins of cos θμfrom fitting T2K CC0π0pdata only. 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 Table III. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-19 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π0pmeasurement 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π0pmeasurement 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π0pNp result) in the forward angle region. Figure 18 additionally shows that the Q2parameter that affects the region 0.25 ≤Q2<0.50 GeV2converges to a value significantly outside of its prior uncertainty when considering the CC0π0pmeasurement. Note that a smaller raising of the parameter was seen the carbon-oxygen measurement considered in Sec. IV B 1. Overall this might indicate a lack of freedom to vary the cross section at highmomentum transfer. The postfit values of the parameters from the simultaneous fit of CC0π0pand CC0πNp 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 δpT 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 δpTmeasurement, we could expect a better separation between the struck nucleon-related uncertainty (FSI, SRC, 2p2h) probed by δpTand the low-energy transfer effects probed mostly by the forward outgoing muon kinematics. 3. Fit to MINERvA cross section data in CC0πNp topology on hydrocarbon Results from the fit to the MINERvA measurement of the CC0πNp (considering Nprotons above 500 MeV=c) cross section as a function of δpTare shown in Fig. 20. FIG. 18. Prefit and postfit values and constraints on the uncertainties from the fit to T2K CC0π0pmeasurement of lepton kinematics and CC0πNp measurement of δpTon 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 Table II. FIG. 17. Prefit (red) and postfit (blue) distributions of δpTfrom fitting the T2K CC0πNp data only. 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 Table III. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-20 As suggested by the chi-square values quoted in Table 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 initial-state nucleon 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 δpTdistribution. 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. 4and 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 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. Figure 22 reports the postfit correlation matrix for this fit. In comparison with the T2K δpTfit, we can notice that the anticorrelation between the p- and the s-shell normalization FIG. 19. Postfit correlation matrices from the fit to T2K CC0π0pmeasurement of lepton kinematics (top left) and CC0πNp measurement of δpT(top right), as well as the simultaneous fit (bottom). PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-21 parameters is less prominent thanks to the finer binning which alleviates their degeneracy with a more precise probe of the shape of the δpTdistribution. Besides, the SRC, 2p2h, FSI, and CCRES uncertainties are correlated as expected since they affect the same high-δpTregion. 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 EQE νbias. Figure 23 shows the prefit and postfit spectra and constraints for the distribution of true neutrino energy and the EQE ν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 EQE ν 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 wellconstrained 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 FIG. 21. Prefit (red) and postfit (blue) values and constraints on the uncertainties from the fit to MINERvA CC0πNp measurement of δpTon 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 Table II. FIG. 22. Postfit correlation matrix from the fit to MINERvA CC0πNp measurement of δpTon carbon. FIG. 20. Prefit (red) and postfit (blue) distributions from the fit to MINERvA CC0πþNp measurement of δpTon 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 Table III. J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-22 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 nbdimensions, where nb¼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 nb-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 measurements show that this parametrization is able to offer a well-motivated means to improve the agreement with FIG. 23. Prefit (red) and postfit (blue) constraints on the true neutrino energy (top) and the bias of EQE ν(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. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-23 respect to the nominal model predictions, especially at the T2K energies. Whilst the fits can directly be used as a model “tune”, their primary purpose within this analysis is to validate the proposed parametrization’s ability to describe cross section measurements. 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 δpT. Within these fits the modification of the SF parameters remains consistent with the nominal SF extracted from electron scattering data, thus implying reasonable agreement with it. However, 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 [55,56], 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. Reference [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 semiinclusive cross section measurements suggest some extensions to the parametrization, in particular related to freedoms within the FSI model, will likely be required. ACKNOWLEDGMENTS S. D. would like to thank Phill Rodrigues for insightful discussions concerning PPP. This work was supported by P2IO LabEx (ANR-10-LABX-0038—Project “BSMNu”)in the framework “Investissements d’Avenir”(ANR-11-IDEX- 0003-01), managed by the Agence Nationale de la Recherche (ANR), France; the Ministry of Science and Higher Education (2023/WK/04) and the National Science Centre (UMO-2018/ 30/E/ST2/00441), Poland; CNRS/IN2P3 and CEA, France; the University of Tokyo ICRR’s Inter-University Research Program FY2023 (Ref. No. 2023i-J-001), Japan; the Spanish Ministerio de Ciencia, Innovación y Universidades and ERDF (European Regional Development Fund) under Contract No. PID2020–114687GB-100 and by the Junta de Andalucía Grant No. FQM160, Spain; the Department of Energy, Office of Science, Office of High Energy Physics, Award No. DE-SC0008475; JSPS KAKENHI Grant No. JP18H05536, USA; DFG research grant (Project No. 517206441), Germany; the Royal Society, Grant No. URF\R1\211661, and the STFC, United Kingdom. [1] L. Alvarez-Ruso et al.,Prog. Part. Nucl. Phys. 100, 1 (2018). [2] K. Abe et al. (T2K Collaboration), Nature (London) 580, 339 (2020);583, E16 (2020). [3] M. A. Acero et al. (NOvA Collaboration), Phys. Rev. D 98, 032012 (2018). [4] K. Abe et al. (Hyper-Kamiokande Collaboration), arXiv: 1805.04163. [5] B. Abi et al. (DUNE Collaboration), J. Instrum. 15, T08008 (2020). [6] K. Abe et al. (T2K Collaboration), arXiv:1901.03750. [7] V. Hewes et al. (DUNE Collaboration), Instruments 5,31 (2021). [8] O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A579, 493 (1994). [9] Y. Hayato and L. Pickering, Eur. Phys. J. Special Topics 230, 4469 (2021). [10] K. Abe et al. (T2K Collaboration), Eur. Phys. J. C 83, 782 (2023). [11] K. Abe et al. (T2K Collaboration), Phys. Rev. D 98, 032003 (2018). [12] K. Abe et al. (T2K Collaboration), Phys. Rev. D 101, 112004 (2020). [13] K. Abe et al. (T2K Collaboration), Phys. Rev. D 108, 112009 (2023). J. CHAKRANI et al. PHYS. REV. D 109, 072006 (2024) 072006-24 [14] K. Abe et al. (T2K Collaboration), Phys. Rev. D 101, 112001 (2020). [15] S. Dolan et al.,Phys. Rev. D 105, 032010 (2022). [16] G. D’Agostini, Nucl. Instrum. Methods Phys. Res., Sect. A 346, 306 (1994). [17] R. W. Peelle, Internal Memorandum, Oak Ridge National Laboratory, Washington DC, USA (unpublished). [18] J. Żmuda, K. M. Graczyk, C. Juszczak, and J. T. Sobczyk, Acta Phys. Pol. B 46, 2329 (2015). [19] C. Andreopoulos et al.,Nucl. Instrum. Methods Phys. Res., Sect. A 614, 87 (2010). [20] D. Dutta et al. (JLab E91013 Collaboration), Phys. Rev. C 68, 064603 (2003). [21] B. Bourguille, J. Nieves, and F. Sánchez, J. High Energy Phys. 04 (2021) 004. [22] S. Dolan, G. D. Megias, and S. Bolognesi, Phys. Rev. D 101, 033003 (2020). [23] J. Schwehr, D. Cherdack, and R. Gran, arXiv:1601.02038. [24] S. Dytman, Y. Hayato, R. Raboanary, J. Sobczyk, J. Tena Vidal, and N. Vololoniaina, Phys. Rev. D 104, 053006 (2021). [25] R. González-Jim´enez, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, N. Jachowicz, G. D. Megias, K. Niewczas, A. Nikolakopoulos, J. W. Van Orden, and J. M. Udías, Phys. Rev. C 105, 025502 (2022). [26] J. M. Franco-Patino, R. González-Jim´enez, S. Dolan, M. B. Barbaro, J. A. Caballero, G. D. Megias, and J. M. Udias, Phys. Rev. D 106, 113005 (2022). [27] A. M. Ankowski, O. Benhar, and M. Sakuda, Phys. Rev. D 91, 033005 (2015). [28] A. Nikolakopoulos, R. González-Jim´enez, N. Jachowicz, K. Niewczas, F. Sánchez, and J. M. Udías, Phys. Rev. C 105, 054603 (2022). [29] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev. C83, 045501 (2011). [30] M. B. Avanzini et al.,Phys. Rev. D 105, 092004 (2022). [31] S. Dolan, arXiv:1810.06043. [32] X. G. Lu, L. Pickering, S. Dolan, G. Barr, D. Coplowe, Y. Uchida, D. Wark, M. O. Wascko, A. Weber, and T. Yuan, Phys. Rev. C 94, 015503 (2016). [33] X. G. Lu et al. (MINERvA Collaboration), Phys. Rev. Lett. 121, 022504 (2018). [34] P. Abratenko et al. (MicroBooNE Collaboration), Phys. Rev. D 108, 053002 (2023). [35] L. Lapikas, Nucl. Phys. A553, 297c (1993). [36] J. M. Franco-Patino, S. Dolan, R. González-Jim´enez, M. B. Barbaro, J. A. Caballero, and G. D. Megias, Phys. Rev. D 109, 013004 (2024). [37] T. Aumann et al.,Prog. Part. Nucl. Phys. 118, 103847 (2021). [38] J. E. Amaro, M. B. Barbaro, J. A. Caballero, R. González- Jim´enez, G. D. Megias, and I. Ruiz Simo, J. Phys. G 47, 124001 (2020). [39] O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R. Seki, Phys. Rev. D 72, 053005 (2005). [40] E. J. Moniz, I. Sick, R. R. Whitney, J. R. Ficenec, R. D. Kephart, and W. P. Trower, Phys. Rev. Lett. 26, 445 (1971). [41] C. Maieron, T. W. Donnelly, and I. Sick, Phys. Rev. C 65, 025502 (2002). [42] A. S. Meyer, A. Walker-Loud, and C. Wilkinson, Annu. Rev. Nucl. Part. Sci. 72, 205 (2022). [43] B. Bhattacharya, R. J. Hill, and G. Paz, Phys. Rev. D 84, 073006 (2011). [44] A. S. Meyer, M. Betancourt, R. Gran, and R. J. Hill, Phys. Rev. D 93, 113015 (2016). [45] K. Niewczas and J. T. Sobczyk, Phys. Rev. C 100, 015505 (2019). [46] A. Ershova et al.,Phys. Rev. D 106, 032009 (2022). [47] J. Isaacson, W. I. Jay, A. Lovato, P. A. N. Machado, and N. Rocco, Phys. Rev. C 103, 015502 (2021). [48] D. Rein and L. M. Sehgal, Ann. Phys. (N.Y.) 133,79 (1981). [49] D. Rein and L. M. Sehgal, Phys. Lett. 104B, 394 (1981); 106B, 513(E) (1981). [50] K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77, 053003 (2008). [51] K. M. Graczyk, J. Żmuda, and J. T. Sobczyk, Phys. Rev. D 90, 093001 (2014). [52] C. Berger and L. M. Sehgal, Phys. Rev. D 76, 113004 (2007). [53] G. M. Radecky et al.,Phys. Rev. D 25, 1161 (1982). [54] T. Kitagaki et al.,Phys. Rev. D 34, 2554 (1986). [55] C. Bronner, Accelerator neutrino I, recent results from T2K (2022), 10.5281/zenodo.6683821. [56] K. Abe et al. (T2K Collaboration), Constraining neutrino interaction uncertainties for neutrino oscillation measurements at the T2K experiment (2023). [57] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B707, 72 (2012). [58] G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. D 94, 013012 (2016). [59] K. Abe et al. (T2K Collaboration), Phys. Rev. D 96, 092006 (2017);98, 019902(E) (2018). [60] K. Abe et al. (T2K Collaboration), Phys. Rev. D 103, 112008 (2021). [61] K. Abe et al. (T2K Collaboration), Phys. Rev. D 101, 112004 (2020). [62] P. Stowell et al.,J. Instrum. 12, P01016 (2017). [63] F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975). [64] P. Abratenko et al. (MicroBooNE Collaboration), Phys. Rev. D 105, 072001 (2022). [65] J. Chakrani, M. Buizza Avanzini, and S. Dolan, Proc. Sci., NuFact2021 (2022) 235 [arXiv:2202.03219]. [66] L. Koch, Phys. Rev. D 103, 113008 (2021). [67] J. Tena-Vidal et al. (GENIE Collaboration), Phys. Rev. D 106, 112001 (2022). [68] T. Katori, A measurement of the muon neutrino charged current quasielastic interaction and a test of Lorentz violation with the MiniBooNE experiment, Ph.D. thesis, Indiana U., 2008. [69] L. Alvarez-Ruso, Proc. Sci., NOW2022 (2023) 026. [70] K. Abe et al. (T2K Collaboration), arXiv:1901.03750. PARAMETRIZED UNCERTAINTIES IN THE SPECTRAL …PHYS. REV. D 109, 072006 (2024) 072006-25