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

J. Chakrani<EMAIL_ADDRESS>Laboratoire Leprince-Ringuet,
CNRS, Ecole polytechnique, Institut Polytechnique de Paris, Palaiseau, France
S. Dolan<EMAIL_ADDRESS>European Organization for Nuclear
Research (CERN), Geneva, Switzerland  M. Buizza Avanzini<EMAIL_ADDRESS>Laboratoire Leprince-Ringuet, CNRS, Ecole polytechnique, Institut
Polytechnique de Paris, Palaiseau, France  A. Ershova  IRFU, CEA, Université
Paris-Saclay, Gif-sur-Yvette, France  L. Koch  Institut für Physik, Johannes
Gutenberg-Universität, Mainz, Germany  K. McFarland  University of Rochester,
Department of Physics and Astronomy, Rochester, New York, U.S.A  G. D. Megias
Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla,
41080 Sevilla, Spain  L. Munteanu  European Organization for Nuclear Research
(CERN), Geneva, Switzerland  L. Pickering  Royal Holloway University of
London, Department of Physics, Egham, Surrey, United Kingdom  K. Skwarczynski
National Centre for Nuclear Research, Warsaw, Poland  V. Q. Nguyen
Laboratoire Leprince-Ringuet, CNRS, Ecole polytechnique, Institut
Polytechnique de Paris, Palaiseau, France  C. Wret  Oxford University, Oxford,
United Kingdom

###### Abstract

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 Alvarez-Ruso  _et al._ (2018).
Whilst statistical uncertainties remain large in the ongoing accelerator-based
long-baseline neutrino oscillation experiments, T2K Abe  _et al._ (2020a) and
NO$\nu$A Acero  _et al._ (2018), the future experiments, Hyper-Kamiokande (HK)
Abe  _et al._ (2018a) and DUNE Abi  _et al._ (2020), are likely to be
dominated by uncertainties related to neutrino interactions. To tackle this
issue, experiments are designing and constructing more sophisticated near
detectors Abe  _et al._ (2019a); Hewes  _et al._ (2021); Abe  _et al._
(2018a). 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_{\nu}\sim 1~{}\textrm{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 Alvarez-Ruso  _et al._ (2018).
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 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. Benhar  _et al._
(1994). 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
Hayato and Pickering (2021), the SF model is used in the T2K experiment’s
recent measurements of neutrino oscillations Abe  _et al._ (2023a) and
neutrino-nucleus cross sections Abe  _et al._ (2018b, 2020b, 2023b, 2020c). It
has also been used to study the sensitivity of T2K’s near-detector upgrade
Dolan  _et al._ (2022).

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) D’Agostini (1994); ppp  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 Hayato and Pickering (2021), NuWro Żmuda  _et al._ (2015) and
GENIE Andreopoulos  _et al._ (2010), 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 “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. Dutta  _et al._ (2003), for CCQE neutrino-nucleus interactions with a
single nucleon in the final state before consideration of FSI (i.e.
$\nu+A\rightarrow\ell+N+A^{\prime}$, where $A$ and $A^{\prime}$ are the
initial-state nucleus and the final-state remnant respectively and $N$ is the
outgoing nucleon), these can be defined as111Whilst removal and missing
energy, and initial state and missing momentum, are equivalent here, it should
be noted that this is for the special case of single nucleon knock out with no
FSI considered within this paper. The missing energy and momentum observables
reported by electron scattering experiments (as in Ref. Dutta  _et al._
(2003)) must usually be corrected in order to obtain removal energy and
initial state momentum.:

$E_{m}=E_{\nu}-E_{\ell}-T_{N}-T_{A^{\prime}}-\Delta m_{N},$ (1)
$p_{m}=|\vec{p}_{m}|=|\vec{p}_{\nu}-\vec{p}_{\ell}-\vec{p}_{N}|,$ (2)

where $E_{\nu}$ and $\vec{p}_{\nu}$ ($E_{\ell}$ and $\vec{p}_{\ell}$) 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^{\prime}}$ is the kinetic energy of the
nuclear remnant; and $\Delta m_{N}$ accounts for the mass difference of the
initial- and final-state nucleon. $T_{A^{\prime}}$ 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}\sim 18$ MeV and a diffuse $s$-shell
at $E_{m}\sim 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
Bourguille  _et al._ (2021)).

Figure 1: Simulations of the two-dimensional distribution of missing energy,
$E_{m}$, and missing momentum, $p_{m}$, for carbon using the RFG (blue), LFG
(green) and SF (red) model simulated with the NEUT neutrino interaction
generator.

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,
$\vec{p}$, below the Fermi momentum $p_{F}$, which yields a cliff feature at
$p_{m}\sim 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.

Figure 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.

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 nucleon222In the case of single-
track events intended to isolate single pion events, the outgoing nucleon mass
can be replaced by the $\Delta(1232)$ mass.,

$E_{\nu}^{\mathrm{QE}}=\frac{2E_{l}\tilde{m}_{N}-\left(m_{l}^{2}+\tilde{m}_{N}^{2}-m_{N}^{2}\right)}{2\left(\tilde{m}_{N}-E_{\ell}+p_{l}\cos\theta_{l}\right)}$
(3)

where $m_{N}$ is the mass of the struck nucleon, $\tilde{m}_{N}=m_{N}-E_{b}$,
and $m_{l}$, $E_{l}$, $p_{l}$, and $\theta_{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_{\nu}$, 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 $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_{\nu}^{\mathrm{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 $\vec{p}_{m}$ of nucleons within the nucleus is
written as:

$P(\vec{p}_{m},E_{m})=P_{\mathrm{MF}}(\vec{p}_{m},E_{m})+P_{\mathrm{corr}}(\vec{p}_{m},E_{m}),$
(4)

where $P_{\mathrm{MF}}(\vec{p}_{m},E_{m})$ is the mean field (MF) part and
$P_{\mathrm{corr}}(\vec{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 $\vec{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 Dolan
_et al._ (2019); Schwehr  _et al._ (2016). 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.

Figure 3: The prediction of the two-dimensional distribution of missing energy
and missing 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 $E_{m},p_{m}$) from the
SRC region (high $E_{m},p_{m}$) in the NuWro (dashed) and NEUT (solid)
neutrino interaction generators.

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 Dytman  _et al._
(2021) and/or a correction to the cross section based on the distortion of the
outgoing nucleon via a nuclear potential (see e.g. González-Jiménez  _et al._
(2022); Franco-Patino  _et al._ (2022)). 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(\vec{p}_{m},E_{m})\rightarrow
P(\vec{p}_{m},E_{m})\theta\left(|\vec{p}_{m}+\vec{q}|-{p}_{F}\right),$ (5)

where $\theta$ is the Heaviside step function, $\vec{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}\lesssim 70$ MeV)—which constitute a significant
fraction ($\sim 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
Dytman  _et al._ (2021) with re-interaction probabilities derived from data.
Another approach to considering the impact of FSI on an SF model is presented
in Ref. Ankowski  _et al._ (2015), 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. Nikolakopoulos  _et al._ (2022). 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 Nieves  _et al._ (2011),
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) Avanzini  _et al._ (2022); Dolan (2018). For
instance, the single-transverse variables Lu  _et al._ (2016) allow a probe of
nuclear effects, particularly with the measurement of the transverse momentum
imbalance $\delta p_{\textrm{T}}$ defined from the kinematics of the outgoing
charged lepton and nucleon as:

$\delta\vec{p}_{\mathrm{T}}=\vec{p}^{\mathrm{T}}_{l}+\vec{p}^{\mathrm{T}}_{N},$
(6)

where $\vec{p}^{\mathrm{T}}_{l}$ and $\vec{p}^{\mathrm{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 $\delta p_{\textrm{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 Abe  _et al._ (2018b) and MINERvA Lu  _et al._ (2018)
measurements of $\delta p_{\textrm{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 Abratenko  _et al._ (2023)
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 Abe  _et al._ (2023a). 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. 2 summarizes all the parametrized systematic uncertainties discussed in
this section, along with their central values and prior uncertainties.

###  III.1 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.

####  III.1.1 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 González-Jiménez  _et al._ (2022); Franco-Patino  _et al._ (2022).
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:

$\displaystyle f_{\mathrm{shell}}(E_{m})$ $\displaystyle=$ $\displaystyle
1+N_{\mathrm{shell}}\times\exp\left(-\frac{(E_{m}-E_{\mathrm{shell}})^{2}}{2\sigma_{\mathrm{shell}}^{2}}\right)$
(7) $\displaystyle=$ $\displaystyle 1+N_{\mathrm{shell}}\times
g_{\mathrm{shell}}(E_{m})$ (8)

where $N_{\mathrm{shell}}$ is the normalization parameter of a given shell,
and $E_{\mathrm{shell}}$ and $\sigma_{\mathrm{shell}}$ correspond to the
center and the width of the Gaussian function, $g_{\mathrm{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_{\mathrm{shell}}$ and $\sigma_{\mathrm{shell}}$
are derived from an analysis of the missing energy distributions in NEUT, and
are shown for each shell in Tab. 1. The impact of varying $N_{\mathrm{shell}}$
for the carbon $p$\- and $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_{\mathrm{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_{\mathrm{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_{\mathrm{shell}}$
to give the event a significant weighting.

Target | Shell |  $E_{\mathrm{shell}}$ [MeV] |  $\sigma_{\mathrm{shell}}$ [MeV] |  $N_{\mathrm{shell}}$ uncertainty
---|---|---|---|---
Carbon | $p$ | 18 | 15 | 0.2
$s$ | 36 | 25 | 0.4
Oxygen | $p_{1/2}$ | 12 | 8 | 0.45
$p_{3/2}$ | 19 | 8 | 0.25
$s$ | 42 | 25 | 0.75
Table 1: 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.

Figure 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_{\mathrm{shell}}$) from
$-1.5\sigma$ (pink) to $+1.5\sigma$ (cyan), compared to the nominal (black).
The interactions were generated on a carbon target using the T2K flux. Figure
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.

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 $\delta
p_{\textrm{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_{\mathrm{shell}}$ parameters conservative. The chosen values are reported
in Tab. 1. These values cover differences in the shell occupancy of the SF
model (derived from $(e,e^{\prime}p)$ data) to those expected in an
independent particle shell model González-Jiménez  _et al._ (2022); Franco-
Patino  _et al._ (2022), whilst ensuring the same variation for each parameter
changes the total CCQE cross section by the same amount ($\sim 10\%$). Fig. 6
shows the impact of varying each shell normalization parameter
$N_{\mathrm{shell}}$ on the total CCQE cross section, and Fig. 4 demonstrates
potential sensitivity to variations through measurements of $\delta
p_{\textrm{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.

Figure 6: The impact of the shell normalization parameters
($N_{\mathrm{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 $\pm 10\%$ variations, chosen to
correspond to the $1\sigma$ uncertainty for these parameters. The interactions
were generated using the T2K flux.

####  III.1.2 Missing-momentum shape

Electron scattering data from Ref. Dutta  _et al._ (2003) 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^{\prime}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]$.

Figure 7: Distributions of the missing momentum from the NEUT SF inputs
(black) compared to electron scattering measurements Dutta  _et al._ (2003)
(blue and red), made for different nucleon and lepton kinematics for carbon in
the $p$-shell (left) and the $s$-shell (right).

These $p_{m}$ shape 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 $\delta p_{\textrm{T}}$ is affected by
this uncertainty since it corresponds, in the absence of other nuclear
effects, to the transverse projection of the Fermi motion.

Figure 8: The distributions for the angle between the outgoing lepton and
neutrino (left) and missing transverse momentum (right) as given by the SF
model implemented in NEUT, showing the impact of variations of the $p$-shell
shape uncertainty from $-2\sigma$ (pink) to $2\sigma$ (cyan), compared to the
nominal (black). The interactions were generated on a carbon target using the
T2K flux.

####  III.1.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. Benhar  _et al._ (1994) 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 $\sim 5\%$ of the total CCQE interactions for both carbon and
oxygen. Alternatively, the SF implementation in NuWro takes a different
approach by making non-rectangular 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
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 $\sim$15% of the total CCQE interactions.

Figure 9: Distribution of the missing momentum in the SF model in NEUT (left)
and NuWro (right) for carbon with the MF and SRC contributions.

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.

###  III.2 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.

####  III.2.1 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
Hayato and Pickering (2021). 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 $\sim$20 MeV/$c$ Benhar  _et al._ (2005); Ankowski  _et
al._ (2015); Moniz  _et al._ (1971); Maieron  _et al._ (2002), 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 $\pm 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.3.2.

Figure 10: The distributions of the angle between the outgoing muon and
neutrino (top) and the muon momentum in the $\cos\theta_{\mu}>0.9$ region
(bottom), showing the impact of varying the PB threshold from $-1.5\sigma$
(cyan) to $+1.5\sigma$ (pink), compared to the nominal (black). The
interactions were generated on a carbon target using the T2K flux.

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 transfers. 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}\sim
450~{}\textrm{MeV}/c$).

####  III.2.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 Nikolakopoulos
_et al._ (2022). To account for this missing alteration to the inclusive cross
section due to FSI–which is physics beyond the PWIA–Ref. Ankowski  _et al._
(2015) 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},\left|\vec{q}\right|)$. 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 $\pm 50\%$ uncertainty. Hence,
the prior acts to stabilize the fit whilst being conservative enough that any
constraint in the fit is dominated by those implied 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: for example, their impact on muon momentum for fixed
ranges of muon angle are different.

Figure 11: The distributions of the angle between the outgoing muon and
neutrino (top) and the muon momentum in the $\cos\theta_{\mu}>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).

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},\left|\vec{q}\right|)$ from NuWro only alters kinematics of the
outgoing nucleon via the lepton-nucleon 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 uncertainties in the
FSI cascade is applied, which allows some freedom to mitigate this issue and
is discussed in Sec. III.3.2.

###  III.3 Additional CC0$\pi$ 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$\pi$, 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_{\nu}\sim
0.6~{}\textrm{GeV}$, and MINERvA has a wide-band beam with average neutrino
energy $\langle E_{\nu}\rangle\sim 3.5~{}\textrm{GeV}$.

####  III.3.1 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_{A}^{\mathrm{QE}}$ in the dipole form factor is set to $1.03\pm
0.06~{}\textrm{GeV}/c^{2}$, estimated from analysis of bubble chamber data Abe
_et al._ (2023a). Additionally, the dipole approximation of the form factor
may be insufficient to describe the evolution of the cross section,
particularly at high $Q^{2}$ Meyer  _et al._ (2022). Therefore, three “high
$Q^{2}$ parameters”, which vary the normalization of CCQE interactions in the
ranges $Q^{2}=[0.25-0.50],[0.50-1.00]$, and $[1.00-\infty]$ in GeV2, are
applied to allow freedom beyond the dipole form factors. This approach is
derived from the $Q^{2}>0.25$ GeV2 parameters in Ref. Abe  _et al._ (2023a).
As also detailed in Ref. Abe  _et al._ (2023a), 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
Bhattacharya  _et al._ (2011); Meyer  _et al._ (2016).

####  III.3.2 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 cross-section measurements are performed in variables that use both
the lepton and the hadron information, such as $\delta p_{\textrm{T}}$. To
account for nucleon FSI uncertainty in this work, CC$0\pi$ 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
Niewczas and Sobczyk (2019). 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 Nikolakopoulos  _et al._ (2022); Ershova  _et al._
(2022); Isaacson  _et al._ (2021)), but is beyond the scope of this work.

####  III.3.3 Resonant pion production

A non-negligible fraction of CC resonant production (CCRES) contributes to the
CC0$\pi$ 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 Rein and Sehgal
(1981a, b); Graczyk and Sobczyk (2008); Graczyk  _et al._ (2014); Berger and
Sehgal (2007) implemented in NEUT. The parameters of this model are: the axial
mass, $M_{A}^{\mathrm{RES}}$, the value of the axial form factor when
$Q^{2}=0$, $C_{5}^{A}$, 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 Radecky  _et al._ (1982) and BNL Kitagaki  _et al._ (1986),
which has been used in Ref. Bronner (2022); t2k .

####  III.3.4 2p2h interactions

As already mentioned, 2p2h interactions also populate the CC0$\pi$ topology.
NEUT describes 2p2h interactions using the Nieves et al. model Nieves  _et
al._ (2012), whose cross section features two distinct peaks in the energy and
momentum transfer that broadly correspond to $\Delta$ and non-$\Delta$
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 Megias  _et al._ (2016); Dolan  _et al._ (2019).
Additionally, we prescribe a 2p2h shape uncertainty which varies the relative
strength between the $\Delta$ and non-$\Delta$ contributions, as used and
described in Refs. Abe  _et al._ (2023a, 2017, 2021).

Parameter | Central value | Prior uncertainty ($1\sigma$) | Notes
---|---|---|---
Carbon parameters
$p$-shell norm. C | 0 | 20% | $1\sigma$ variation changes CCQE cross section by 10%
$s$-shell norm. C | 0 | 40%
$p$-shell shape C | 0 | 100% | From $(e,e^{\prime}p)$ data, used only for fits in $\delta p_{\textrm{T}}$ (CH target)
$s$-shell shape C | 0 | 100%
SRC norm. C | 1 | 100% |
Pauli Blocking C (MeV/$c$) | 209 | 30 | Used only with fits in lepton kinematics
Optical Potential C | 50% | 50%
2p2h norm. C | 1 | 30% |
2p2h shape C | 0 | 300% | Defined in Ref. Abe  _et al._ (2021)
Pion abs. norm. C | 1 | 30% |
Oxygen parameters
$p_{1/2}$-shell norm. O | 0 | 45% | $1\sigma$ variation changes CCQE cross section by 10%
$p_{3/2}$-shell norm. O | 0 | 25%
$s$-shell norm. O | 0 | 75%
SRC norm. O | 1 | 100% |
Pauli Blocking O (MeV/$c$) | 209 | 30 | Used only with fits in lepton kinematics
Optical Potential O | 50% | 50%
2p2h norm. O | 1 | 30% |
2p2h shape O | 0 | 300% | Defined in Ref. Abe  _et al._ (2021)
Pion abs. norm. O | 1 | 30% |
Nucleon interaction parameters
$M_{A}^{\mathrm{QE}}$ (GeV/$c^{2}$) | 1.03 | 0.06 | Dipole parametrization
High $Q^{2}$ norm. 1 | 1 | 11% | $Q^{2}\in[0.25,0.50\text{ GeV}^{2}[$
High $Q^{2}$ norm. 2 | 1 | 18% | $Q^{2}\in[0.50,1.00\text{ GeV}^{2}[$
High $Q^{2}$ norm. 3 | 1 | 40% | $Q^{2}\in[1.00\text{ GeV}^{2},+\infty[$
$M_{A}^{\mathrm{RES}}$ (GeV/$c^{2}$) | 0.91 | 0.1 | Correlated Bronner (2022); t2k
$C_{5}^{A}(Q^{2}=0)$ | 1.06 | $10\%$
$I_{1/2}$ non-res. bkg. | 1.21 | $27\%$
FSI parameters
With nucleon FSI | 1 | 30% |  Strongly anti-correlated, see Sec. III.3.2 Used only in fits to semi-inclusive measurements333I.e. measurements with restrictions on the outgoing hadron kinematics, corresponding to the three fits of Sec. IV.2.2 and the fit of Sec. IV.2.3
Without nucleon FSI | 1 | 30%
Table 2: Summary of the parameters introduced in Sec. III and their prior
uncertainties.

##  IV Fits to CC$0\pi$ 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$\pi$, 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 cross-section measurements
reported in muon kinematics and the transverse kinematic imbalance between the
outgoing lepton and nucleons, such as $\delta p_{\textrm{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$\pi$ cross section in muon kinematics
$(p_{\mu},\cos\theta_{\mu})$ on oxygen and carbon Abe  _et al._ (2020d), the
T2K CC0$\pi$ cross section in $\delta p_{\textrm{T}}$, muon, and proton
kinematics on CH Abe  _et al._ (2018b), and the MINERvA CC0$\pi$ cross section
in $\delta p_{\textrm{T}}$ on CH Lu  _et al._ (2018). 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. 2 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 Stowell  _et al._ (2017). 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 James and
Roos (1975).

###  IV.1 Chi-squared test-statistic

To perform the desired fits, we developed and implemented in NUISANCE a
specific $\chi^{2}$ test-statistic, with the explicit goal of mitigating the
impact of Peelle’s Pertinent Puzzle (PPP) D’Agostini (1994); ppp , 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=\mathrm{Cov}\left[\\{B_{i}\\}\right]$ as provided by the experiments and
minimizing a standard test-statistic like the chi-squared:

$\chi_{\mathrm{}}^{2}=\sum_{i,j}\left(B_{i}-B_{i}^{\mathrm{MC}}\right)\left(M^{-1}\right)_{i,j}\left(B_{j}-B_{j}^{\mathrm{MC}}\right)$
(9)

where $B^{\mathrm{MC}}=(B^{\mathrm{MC}}_{1},...,B^{\mathrm{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 $\chi^{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
$\chi_{\mathrm{}}^{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
D’Agostini (1994). In the case of non-Gaussian uncertainties the $\chi^{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 Abe  _et al._ (2020b); Dolan (2018), or neglecting the bin-to-bin
correlations Abratenko  _et al._ (2022). 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 uncertainty that remains constant as a
function of the normalization, as motivated by the arguments of Ref.
D’Agostini (1994). 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:

$C_{i}=f(B_{i})=\begin{cases}\alpha\frac{B_{i}}{\sum_{k}B_{k}},&1\leq i\leq
n-1\\\ B_{T}=\sum_{k}B_{k},&i=n\end{cases}$ (10)

where $\alpha$ 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\mapsto C$ is bijective, meaning that no information is lost by
moving from $B$ to $C$.

In this new basis, a covariance matrix,
$N=\mathrm{Cov}\left[\\{C_{i}\\}\right]$, 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:

$N=J(f).M.J(f)^{T}$ (11)

where $J(f)$ is the Jacobian of the non-linear transformation $f$. The new
covariance matrix is expressed as follows:

$N=\begin{pmatrix}(N_{\mathrm{S}})_{i,j}=\frac{\alpha^{2}}{B_{T}^{2}}\left(M_{i,j}-\frac{B_{i}}{B_{T}}\sum_{l}M_{i,l}-\frac{B_{j}}{B_{T}}\sum_{k}M_{k,j}-\frac{B_{i}B_{j}}{B_{T}^{2}}\sum_{kl}M_{k,l}\right)&\vline&\begin{matrix}\frac{\alpha}{B_{T}}\left(\sum_{l}M_{1,l}-\frac{B_{1}}{B_{T}}\sum_{kl}M_{k,l}\right)\\\\[10.00002pt]
\vdots\\\\[10.00002pt]
\frac{\alpha}{B_{T}}\left(\sum_{l}M_{n-1,l}-\frac{B_{n-1}}{B_{T}}\sum_{kl}M_{k,l}\right)\end{matrix}\\\\[10.00002pt]
\hline\cr\begin{matrix}\frac{\alpha}{B_{T}}\left(\sum_{k}M_{k,1}-\frac{B_{1}}{B_{T}}\sum_{kl}M_{k,l}\right)&\cdots&\frac{\alpha}{B_{T}}\left(\sum_{k}M_{k,n-1}-\frac{B_{n-1}}{B_{T}}\sum_{kl}M_{k,l}\right)\end{matrix}&\vline&\sum_{kl}M_{k,l}\end{pmatrix}$
(12)

The matrix $N$ has the same dimension and the same positive-definiteness
properties as $M$ since the mapping $B\mapsto C$ is a bijection. $N$ is
composed of two diagonal blocks: the $N_{\mathrm{S}}$ block which corresponds
to the shape-only covariance, and the $\sum_{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:

$\chi_{\mathrm{NS}}^{2}=\sum_{1\leq i,j\leq
n}\left(C_{i}-C_{i}^{\mathrm{MC}}\right)\left(N^{-1}\right)_{i,j}\left(C_{j}-C_{j}^{\mathrm{MC}}\right).$
(13)

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 $\chi_{\mathrm{NS}}^{2}$ allows to mitigate the PPP problem observed
when using the standard $\chi^{2}$ of Eq. 9 in the fit. This method was first
discussed in Ref. Chakrani  _et al._ (2022); Koch (2021) and also used, via
our implementation, in Refs. Tena-Vidal  _et al._ (2022); Abratenko  _et al._
(2022)444On a historical note: the MiniBooNE collaboration had developed (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. Katori (2008)..

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 assumptions555If 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.. As such, the modified
test-statistic is used for fits but the usual $\chi^{2}$ from Eq. 9 is also
reported.

###  IV.2 Fit results

In this section, the results of the fits, using the $\chi_{\mathrm{NS}}^{2}$
defined in the previous section, to T2K and MINERvA data are presented and
discussed.

####  IV.2.1 Fit to T2K CC0$\pi$ cross section data on oxygen and carbon in
muon kinematics

Figure 12: Prefit (red) and postfit (blue) distributions of $p_{\mu}$ in bins
of $\cos\theta_{\mu}$ from the fit to T2K CC0$\pi$ 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 $\chi^{2}_{\mathrm{NS}}$
used in the minimization is reported in Tab. 3. Figure 13: Prefit (red) and
postfit (blue) values and constraints on the uncertainties from the fit to T2K
CC0$\pi$ 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. 2.  Figure 14: Postfit correlation matrix from the fit to T2K CC0$\pi$
joint measurement of lepton kinematics on carbon and oxygen.

Fig. 12 shows the prefit and postfit cross-section predictions for the T2K
CC0$\pi$ measurement on carbon and oxygen targets from Ref. Abe  _et al._
(2020d), 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
$\chi_{\mathrm{NS}}^{2}$ as well as the usual $\chi^{2}$ are reported in Tab.
3.

Measurement | Prefit $\chi^{2}_{\mathrm{NS}}$ | Postfit $\chi^{2}_{\mathrm{NS}}$ | Prefit $\chi^{2}$ | Postfit $\chi^{2}$ | Number of bins
---|---|---|---|---|---
T2K oxygen + carbon (Sec. IV.2.1) | 110.88 | 35.81 | 98.79 | 30.30 | 58
T2K CC$0\pi Np$ $\delta p_{\textrm{T}}$ only (Sec. IV.2.2) | 12.59 | 7.37 | 15.72 | 8.48 | 8
T2K CC$0\pi 0p$ $(p_{\mu},\cos\theta_{\mu})$ only (Sec. IV.2.2) | 144.35 | 75.13 | 107.57 | 62.55 | 50
T2K CC$0\pi 0p$ \+ CC$0\pi Np$ (Sec. IV.2.2) | $144.35+14.56$ | $86.80+10.01$ | $107.57+16.76$ | $64.19+11.83$ | $50+8$
MINERvA $\delta p_{\textrm{T}}$ (Sec. IV.2.3) | 109.10 | 79.51 | 114.32 | 76.14 | 24
Table 3: Summary of the prefit and postfit norm-shape chi-square
$\chi^{2}_{\mathrm{NS}}$ used in the minimization, and the usual chi-square
$\chi^{2}$ for reference in the different fits presented in Sec. IV.2 along
with the corresponding number of bins.

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 $\chi_{\mathrm{NS}}^{2}$ and
$\chi^{2}$. Importantly, both are also lower than the number of bins, implying
a quantitatively good agreement between the post-fit model and the measured
cross section. It should be noticed that, as detailed in Ref. Abe  _et al._
(2020d), the prefit data/model disagreement comes mainly from the bins that
correspond to forward lepton kinematics (i.e. the most forward
cos$\theta_{\mu}$ 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 Nieves  _et al._ (2011). 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.2: 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_{A}^{\mathrm{QE}}$ 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
$\sim$250 MeV/$c$. Interpreted as a Fermi-gas Fermi momentum, this would imply
nonphysical nuclear density Alvarez-Ruso (2023) 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_{A}^{\mathrm{QE}}$,
$M_{A}^{\mathrm{RES}}$, 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.

####  IV.2.2 Fit to T2K cross section data in CC0$\pi$0p and CC0$\pi N$p
topologies on hydrocarbon

In this section, we fit T2K CC$0\pi$ cross-section measurements on hydrocarbon
as a function of transverse momentum imbalance $\delta p_{\textrm{T}}$ and
final-state muon kinematics for CC$0\pi$ topologies with (CC$0\pi Np$, $N\geq
1$) and without (CC$0\pi 0p$) a measured proton in the final state
respectively Abe  _et al._ (2018b). Note that $0p$ and $Np$ refer to protons
above and below tracking threshold ($\sim$500 MeV/c). Three distinct fits are
performed:

  1. (a)

fitting lepton kinematics in CC$0\pi 0p$,

  2. (b)

fitting $\delta p_{\textrm{T}}$ in CC$0\pi Np$,

  3. (c)

simultaneously fitting lepton kinematics in CC$0\pi 0p$ and $\delta
p_{\textrm{T}}$ in CC$0\pi 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
CC$0\pi 0p$ and the CC$0\pi Np$ topologies correspond to distinct regions of
momentum transfer (as a proton requires $\sim$450 MeV momentum to be observed)
and therefore to different regions of $E_{\nu}$, 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 CC$0\pi 0p$ topology allows protons up to 500 MeV/$c$,
whilst the $\delta p_{\textrm{T}}$ measurement considers only protons above
450 MeV/$c$ although a more complete analysis would provide correlations from
systematic uncertainties.

Figure 15: True squared four momentum transfer (upper) and neutrino energy
(lower) distributions of interactions falling into the the CC$0\pi 0p$ (blue)
and CC$0\pi Np$ (red) topologies in the T2K beam on a carbon target as
predicted by NEUT.

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 CC$0\pi
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
reconstruction666Whilst these bins are in the data release, they were not
shown within Ref. Abe  _et al._ (2018b)..

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. 2, it should first be noted that the three fits do
not use the same set of parameters.

Figure 16: Prefit (red) and postfit (blue) distributions of $p_{\mu}$ in bins
of $\cos\theta_{\mu}$ from fitting T2K CC0$\pi 0p$ data only. The usual chi-
squares as well as the number of bins are quoted in the legend. The NS chi-
square $\chi^{2}_{\mathrm{NS}}$ used in the minimization is reported in Tab.
3. Figure 17: Prefit (red) and postfit (blue) distributions of $\delta
p_{\textrm{T}}$ from fitting the T2K CC0$\pi Np$ data only. The usual chi-
squares as well as the number of bins are quoted in the legend. The NS chi-
square $\chi^{2}_{\mathrm{NS}}$ used in the minimization is reported in Tab.
3. Figure 18: Prefit and postfit values and constraints on the uncertainties
from the fit to T2K CC0$\pi 0p$ measurement of lepton kinematics and CC0$\pi
Np$ measurement of $\delta p_{\textrm{T}}$ 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. 2.

Figure 19: Postfit correlation matrices from the fit to T2K CC0$\pi 0p$
measurement of lepton kinematics (top left) and CC0$\pi Np$ measurement of
$\delta p_{\textrm{T}}$ (top right), as well as the simultaneous fit (bottom).

The binning of $\delta p_{\textrm{T}}$, as shown in Fig. 17, is rather coarse,
and so even the prefit chi-squares reported in Tab. 3 for the CC$0\pi 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 $\delta
p_{\textrm{T}}$, particularly in the tail of its distribution Lu  _et al._
(2016). 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 $\delta p_{\textrm{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 $\delta p_{\textrm{T}}$ as a function of
$\delta\alpha_{\textrm{T}}$ (another transverse variable Lu  _et al._ (2016)),
or even as a function of outgoing lepton kinematics. Such a measurement has
recently been performed by the MicroBooNE collaboration Abratenko  _et al._
(2023), and there are prospects for similar high-precision measurements using
the upgraded T2K ND280 detector Abe  _et al._ (2019b). On the other hand, even
with the coarse binning, the CCQE dominance in the bulk of $\delta
p_{\textrm{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 CC$0\pi Np$ fit and the oxygen + carbon analysis from Sec.
IV.2.1, fits including the CC$0\pi 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 $\delta p_{\textrm{T}}$, but is insufficient to fully
describe the CC$0\pi 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.2.1, the improved agreement of the model with the data in the CC$0\pi 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 CC$0\pi 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 CC$0\pi 0pNp$ result) in the
forward angle region.

Fig. 18 additionally shows that the $Q^{2}$ parameter that affects the region
$0.25\leq Q^{2}<0.50$ GeV2 converges to a value significantly outside of its
prior uncertainty when considering the CC$0\pi 0p$ measurement. Note that a
smaller raising of the parameter was seen the carbon-oxygen measurement
considered in Sec. IV.2.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 CC$0\pi 0p$
and CC$0\pi Np$ measurements are also displayed in Fig. 18 (purple). Most of
the parameters converge to similar values as in the CC$0\pi 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 CC$0\pi 0p$-only one thanks to the constraint from the
$\delta p_{\textrm{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 CC$0\pi 0p$-only and the simultaneous
fit. With a more statistically significant $\delta p_{\textrm{T}}$
measurement, we could expect a better separation between the struck nucleon-
related uncertainty (FSI, SRC, 2p2h) probed by $\delta p_{\textrm{T}}$ and the
low energy transfer effects probed mostly by the forward outgoing muon
kinematics.

####  IV.2.3 Fit to MINERvA cross section data in CC0$\pi N$p topology on
hydrocarbon

Results from the fit to the MINERvA measurement of the CC$0\pi Np$
(considering $N$ protons above 500 MeV/$c$) cross section as a function of
$\delta p_{\textrm{T}}$ are shown in Fig. 20. As suggested by the chi-square
values quoted in Tab. 3, 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 $\delta
p_{\textrm{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
$\chi^{2}_{\mathrm{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 $\sim
20\%$ of the CC$0\pi$ 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 20: Prefit (red) and postfit (blue) distributions from the fit to
MINERvA CC0$\pi+Np$ measurement of $\delta p_{\textrm{T}}$ on carbon. The
usual chi-squares as well as the number of bins are quoted in the legend. The
NS chi-square $\chi^{2}_{\mathrm{NS}}$ used in the minimization is reported in
Tab. 3. Figure 21: Prefit (red) and postfit (blue) values and constraints on
the uncertainties from the fit to MINERvA CC0$\pi Np$ measurement of $\delta
p_{\textrm{T}}$ 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. 2. Figure 22: Postfit correlation matrix
from the fit to MINERvA CC0$\pi Np$ measurement of $\delta p_{\textrm{T}}$ on
carbon.

Fig. 22 reports the postfit correlation matrix for this fit. In comparison
with the T2K $\delta p_{\textrm{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 $\delta
p_{\textrm{T}}$ distribution. Besides, the SRC, 2p2h, FSI and CCRES
uncertainties are correlated as expected since they affect the same
high-$\delta p_{\textrm{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_{\nu}^{\mathrm{QE}}$ bias.

Fig. 23 shows the prefit and postfit spectra and constraints for the
distribution of true neutrino energy and the $E_{\nu}^{\mathrm{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.2.1.

Figure 23: Prefit (red) and postfit (blue) constraints on the true neutrino
energy (top) and the bias of $E_{\nu}^{\mathrm{QE}}$ (bottom) differential
cross sections in the T2K beam on a carbon target following the fit to the T2K
CC0$\pi$ joint measurement of lepton kinematics on carbon and oxygen (Sec.
IV.2.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.

The postfit constraints on the cross section as a function of $E_{\nu}$ 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_{\nu}^{\mathrm{QE}}$ particularly around the true $E_{\nu}$ (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_{\nu}$ 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$\pi$ cross-section measurements 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 $\delta p_{\textrm{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 Bronner (2022); t2k
, 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 Abe  _et al._ (2019a). 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. Dolan  _et al._ (2022) 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.

## Acknowledgements

SD 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
Education and Science (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. 2023i-J-001),
Japan; the Spanish Ministerio de Ciencia, Innovación y Universidades and ERDF
(European Regional Development Fund) under contract 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 DE-SC0008475; JSPS
KAKENHI Grant Number JP18H05536, USA; DFG research grant (project number
517206441), Germany; the Royal Society, grant number URF\R1\211661, and the
STFC, UK.

## References

  * Alvarez-Ruso  _et al._ (2018) L. Alvarez-Ruso _et al._ , Prog. Part. Nucl. Phys. 100, 1 (2018), arXiv:1706.03621 [hep-ph] .
  * Abe  _et al._ (2020a) K. Abe _et al._ (T2K), Nature 580, 339 (2020a), [Erratum: Nature 583, E16 (2020)], arXiv:1910.03887 [hep-ex] .
  * Acero  _et al._ (2018) M. A. Acero _et al._ (NOvA), Phys. Rev. D98, 032012 (2018), arXiv:1806.00096 [hep-ex] .
  * Abe  _et al._ (2018a) K. Abe _et al._ (Hyper-Kamiokande), (2018a), arXiv:1805.04163 [physics.ins-det] .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE collaboration), JINST 15, T08008 (2020), arXiv:2002.02967 [physics.ins-det] .
  * Abe  _et al._ (2019a) K. Abe _et al._ (T2K), (2019a), arXiv:1901.03750 [physics.ins-det] .
  * Hewes  _et al._ (2021) V. Hewes _et al._ (DUNE), Instruments 5, 31 (2021), arXiv:2103.13910 [physics.ins-det] .
  * Benhar  _et al._ (1994) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A579, 493 (1994).
  * Hayato and Pickering (2021) Y. Hayato and L. Pickering,  (2021), 10.1140/epjs/s11734-021-00287-7, arXiv:2106.15809 [hep-ph] .
  * Abe  _et al._ (2023a) K. Abe _et al._ (T2K), (2023a), arXiv:2303.03222 [hep-ex] .
  * Abe  _et al._ (2018b) K. Abe _et al._ (T2K), Phys. Rev. D 98, 032003 (2018b), arXiv:1802.05078 [hep-ex] .
  * Abe  _et al._ (2020b) K. Abe _et al._ (T2K), Phys. Rev. D 101, 112004 (2020b), arXiv:2004.05434 [hep-ex] .
  * Abe  _et al._ (2023b) K. Abe _et al._ (T2K), (2023b), arXiv:2303.14228 [hep-ex] .
  * Abe  _et al._ (2020c) K. Abe _et al._ (T2K), Phys. Rev. D 101, 112001 (2020c), arXiv:2002.09323 [hep-ex] .
  * Dolan  _et al._ (2022) S. Dolan _et al._ , Phys. Rev. D 105, 032010 (2022), arXiv:2108.11779 [hep-ex] .
  * D’Agostini (1994) G. D’Agostini, Nucl. Instrum. Meth. A 346, 306 (1994).
  * (17) Internal Memorandum, Oak Ridge National Laboratory, Washington DC, USA, unpublished R.W. Peelle.
  * Żmuda  _et al._ (2015) J. Żmuda, K. M. Graczyk, C. Juszczak, and J. T. Sobczyk, Acta Phys. Polon. B 46, 2329 (2015), arXiv:1510.03268 [hep-ph] .
  * Andreopoulos  _et al._ (2010) C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Dutta  _et al._ (2003) D. Dutta _et al._ (JLab E91013), Phys. Rev. C 68, 064603 (2003), arXiv:nucl-ex/0303011 .
  * Bourguille  _et al._ (2021) B. Bourguille, J. Nieves, and F. Sánchez, JHEP 04, 004 (2021), arXiv:2012.12653 [hep-ph] .
  * Dolan  _et al._ (2019) S. Dolan, G. D. Megias, and S. Bolognesi, (2019), arXiv:1905.08556 [hep-ex] .
  * Schwehr  _et al._ (2016) J. Schwehr, D. Cherdack, and R. Gran, (2016), arXiv:1601.02038 [hep-ph] .
  * Dytman  _et al._ (2021) S. Dytman, Y. Hayato, R. Raboanary, J. Sobczyk, J. Tena Vidal, and N. Vololoniaina, (2021), arXiv:2103.07535 [hep-ph] .
  * González-Jiménez  _et al._ (2022) R. González-Jiménez, 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), arXiv:2104.01701 [nucl-th] .
  * Franco-Patino  _et al._ (2022) J. M. Franco-Patino, R. González-Jiménez, S. Dolan, M. B. Barbaro, J. A. Caballero, G. D. Megias, and J. M. Udias, Phys. Rev. D 106, 113005 (2022), arXiv:2207.02086 [nucl-th] .
  * Ankowski  _et al._ (2015) A. M. Ankowski, O. Benhar, and M. Sakuda, Phys. Rev. D91, 033005 (2015), arXiv:1404.5687 [nucl-th] .
  * Nikolakopoulos  _et al._ (2022) A. Nikolakopoulos, R. González-Jiménez, N. Jachowicz, K. Niewczas, F. Sánchez, and J. M. Udías, Phys. Rev. C 105, 054603 (2022), arXiv:2202.01689 [nucl-th] .
  * Nieves  _et al._ (2011) J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Rev. C 83, 045501 (2011), arXiv:1102.2777 [hep-ph] .
  * Avanzini  _et al._ (2022) M. B. Avanzini _et al._ , Phys. Rev. D 105, 092004 (2022), arXiv:2112.09194 [hep-ex] .
  * Dolan (2018) S. Dolan, (2018), arXiv:1810.06043 [hep-ex] .
  * Lu  _et al._ (2016) 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. C94, 015503 (2016), arXiv:1512.05748 [nucl-th] .
  * Lu  _et al._ (2018) X. G. Lu _et al._ (MINERvA), Phys. Rev. Lett. 121, 022504 (2018), arXiv:1805.05486 [hep-ex] .
  * Abratenko  _et al._ (2023) P. Abratenko _et al._ (MicroBooNE), (2023), arXiv:2301.03700 [hep-ex] .
  * Benhar  _et al._ (2005) O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R. Seki, Phys. Rev. D 72, 053005 (2005), arXiv:hep-ph/0506116 .
  * Moniz  _et al._ (1971) E. J. Moniz, I. Sick, R. R. Whitney, J. R. Ficenec, R. D. Kephart, and W. P. Trower, Phys. Rev. Lett. 26, 445 (1971).
  * Maieron  _et al._ (2002) C. Maieron, T. W. Donnelly, and I. Sick, Phys. Rev. C 65, 025502 (2002).
  * Meyer  _et al._ (2022) A. S. Meyer, A. Walker-Loud, and C. Wilkinson,  (2022), 10.1146/annurev-nucl-010622-120608, arXiv:2201.01839 [hep-lat] .
  * Bhattacharya  _et al._ (2011) B. Bhattacharya, R. J. Hill, and G. Paz, Phys. Rev. D 84, 073006 (2011), arXiv:1108.0423 [hep-ph] .
  * Meyer  _et al._ (2016) A. S. Meyer, M. Betancourt, R. Gran, and R. J. Hill, Phys. Rev. D 93, 113015 (2016), arXiv:1603.03048 [hep-ph] .
  * Niewczas and Sobczyk (2019) K. Niewczas and J. T. Sobczyk, Phys. Rev. C 100, 015505 (2019), arXiv:1902.05618 [hep-ex] .
  * Ershova  _et al._ (2022) A. Ershova _et al._ , Phys. Rev. D 106, 032009 (2022), arXiv:2202.10402 [hep-ph] .
  * Isaacson  _et al._ (2021) J. Isaacson, W. I. Jay, A. Lovato, P. A. N. Machado, and N. Rocco, Phys. Rev. C 103, 015502 (2021), arXiv:2007.15570 [hep-ph] .
  * Rein and Sehgal (1981a) D. Rein and L. M. Sehgal, Annals Phys. 133, 79 (1981a).
  * Rein and Sehgal (1981b) D. Rein and L. M. Sehgal, Phys. Lett. 104B, 394 (1981b), [Erratum: Phys. Lett.106B,513(1981)].
  * Graczyk and Sobczyk (2008) K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77, 053003 (2008), arXiv:0709.4634 [hep-ph] .
  * Graczyk  _et al._ (2014) K. M. Graczyk, J. Żmuda, and J. T. Sobczyk, Phys. Rev. D 90, 093001 (2014), arXiv:1407.5445 [hep-ph] .
  * Berger and Sehgal (2007) C. Berger and L. M. Sehgal, Phys. Rev. D 76, 113004 (2007), arXiv:0709.4378 [hep-ph] .
  * Radecky  _et al._ (1982) G. M. Radecky _et al._ , Phys. Rev. D 25, 1161 (1982).
  * Kitagaki  _et al._ (1986) T. Kitagaki _et al._ , Phys. Rev. D 34, 2554 (1986).
  * Bronner (2022) C. Bronner, “Accelerator neutrino i, recent results from t2k,”  (2022).
  * (52) T2K collaboration, paper in preparation.
  * Nieves  _et al._ (2012) J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B707, 72 (2012), arXiv:1106.5374 [hep-ph] .
  * Megias  _et al._ (2016) G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. D 94, 013012 (2016), arXiv:1603.08396 [nucl-th] .
  * Abe  _et al._ (2017) K. Abe _et al._ (T2K), Phys. Rev. D 96, 092006 (2017), [Erratum: Phys.Rev.D 98, 019902 (2018)], arXiv:1707.01048 [hep-ex] .
  * Abe  _et al._ (2021) K. Abe _et al._ (T2K), Phys. Rev. D 103, 112008 (2021), arXiv:2101.03779 [hep-ex] .
  * Abe  _et al._ (2020d) K. Abe _et al._ (T2K), Phys. Rev. D 101, 112004 (2020d), arXiv:2004.05434 [hep-ex] .
  * Stowell  _et al._ (2017) P. Stowell _et al._ , JINST 12, P01016 (2017), arXiv:1612.07393 [hep-ex] .
  * James and Roos (1975) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
  * Abratenko  _et al._ (2022) P. Abratenko _et al._ (MicroBooNE), Phys. Rev. D 105, 072001 (2022), arXiv:2110.14028 [hep-ex] .
  * Chakrani  _et al._ (2022) J. Chakrani, M. Buizza Avanzini, and S. Dolan, PoS NuFact2021, 235 (2022), arXiv:2202.03219 [hep-ph] .
  * Koch (2021) L. Koch, Phys. Rev. D 103, 113008 (2021), arXiv:2102.06172 [physics.data-an] .
  * Tena-Vidal  _et al._ (2022) J. Tena-Vidal _et al._ (GENIE), Phys. Rev. D 106, 112001 (2022), arXiv:2206.11050 [hep-ph] .
  * Katori (2008) 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).
  * Alvarez-Ruso (2023) L. Alvarez-Ruso, PoS NOW2022, 026 (2023).
  * Abe  _et al._ (2019b) K. Abe _et al._ (T2K), (2019b), arXiv:1901.03750 [physics.ins-det] .