# Form factor and model dependence in neutrino-nucleus cross section
predictions

Daniel Simons  Department of Physics and Astronomy, University of Iowa, Iowa
City, IA 52242, USA  Noah Steinberg  Theoretical Physics Department, Fermi
National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60410, USA
Alessandro Lovato  Physics Division, Argonne National Laboratory, Argonne,
Illinois 60439, USA  INFN-TIFPA Trento Institute of Fundamental Physics and
Applications, Via Sommarive, 14, 38123 Trento, Italy  Yannick Meurice
Department of Physics and Astronomy, University of Iowa, Iowa City, IA 52242,
USA  Noemi Rocco  Theoretical Physics Department, Fermi National Accelerator
Laboratory, P.O. Box 500, Batavia, IL 60410, USA  Michael Wagman  Theoretical
Physics Department, Fermi National Accelerator Laboratory, P.O. Box 500,
Batavia, IL 60410, USA

###### Abstract

To achieve its design goals, the next generation of neutrino-oscillation
accelerator experiments requires percent-level predictions of neutrino-nucleus
cross sections supplemented by robust estimates of the theoretical
uncertainties involved. The latter arise from both approximations in solving
the nuclear many-body problem and in the determination of the single- and few-
nucleon quantities taken as input by many-body methods. To quantify both types
of uncertainty, we compute flux-averaged double-differential cross sections
using the Green’s function Monte Carlo and spectral function methods as well
as different parameterizations of the nucleon axial form factors based on
either deuterium bubble-chamber data or lattice quantum chromodynamics
calculations. The cross-section results are compared with available
experimental data from the MiniBooNE and T2K collaborations. We also discuss
the uncertainties associated with $N\rightarrow\Delta$ transition form factors
that enter the two-body current operator. We quantify the relations between
neutrino-nucleus cross section and nucleon form factor uncertainties. These
relations enable us to determine the form factor precision targets required to
achieve a given cross-section precision.

††preprint: FERMILAB-PUB-22-745-T

##  I Introduction

The study of neutrino processes is driven by deep questions whose answers may
profoundly change our understanding of physics. In particular, given that
neutrinos have mass and mix, the accelerator-neutrino program aims at
precisely measuring the parameters that characterize their oscillations,
investigating the possible existence of a fourth neutrino flavor, and testing
additional Beyond the Standard Model scenarios. The success of these
experiments rests on our ability to compute neutrino-nucleus cross sections
with quantified theoretical uncertainties [1]. The latter presently yield a
sizable contribution to the total error budget of oscillation parameters [2,
3].

On the other hand, accelerator neutrino experiments allow us access to aspects
of nuclear dynamics that would otherwise be difficult to probe at electron-
scattering facilities. The chief example is the axial form factor of the
nucleon. Its experimental determination dates back to bubble chamber
experiments carried out in the 70s and 80s [4, 5, 6, 7, 8] and to electroweak
single pion production measured at ANL and BNL in the 80s [9, 10, 11, 12]. A
simple dipole parameterization with axial mass $M_{A}\sim 1$ GeV reproduces
single-nucleon data. In contrast, a larger value $M_{A}\sim 1.2$ GeV was
required in order to make relativistic Fermi Gas predictions for the
neutrino-12C cross sections compatible with MiniBooNE data [13]. Note that the
spectral-function formalism, which includes the vast majority of nuclear
correlations, requires an even larger value of $M_{A}$ to reproduce
experimental data [14].

This apparent inconsistency appeared to be solved by models that include two-
body current operators, as they can reproduce MiniBooNE and T2K data with
$M_{A}~{}\sim~{}1$ GeV [15, 16, 17, 18, 19, 20, 21]. However, these models are
somewhat simplified, as they are based on a mean-field description of nuclear
dynamics. As in the one-body current case, it may well be that when nuclear
correlations are accounted for there is room for larger values of $M_{A}$.
This is confirmed by recent Green’s function Monte Carlo (GFMC) results, which
provide a full account of nuclear correlations and two-body current effects
[22]. However, the non-relativistic nature of the GFMC hampers its
applicability to neutrino accelerator experiments where the neutrino flux
energy is of the order of a few GeV. The spectral function (SF) method,
relying on the factorization of the final hadronic state, allows for the
inclusion of relativistic effects and exclusive channels while retaining most
of the important effects coming from multi-nucleon dynamics. The factorization
scheme has been extended to include one- and two-body current operators in a
consistent fashion as well as pion production amplitudes and validated against
electron scattering data [23, 24]. The spectral function of light and medium
mass nuclei has been recently computed exploiting quantum Monte Carlo (QMC)
techniques and shares with the GFMC the same description of nuclear dynamics
[25, 26]. Comparing the results obtained for lepton-nucleus scattering using
these two different approaches enables a precise quantification of the
uncertainties inherent to factorization schemes that need to be accounted for
when assessing the total error of the theoretical calculations in neutrino
oscillation analysis.

There has been significant recent progress in lattice quantum chromodynamics
(LQCD) calculations of nucleon axial and vector form factors [27, 28, 29, 30,
31, 32, 33, 34, 35], which have now been performed using approximately
physical values of the quark masses as well as multiple lattice spacings and
volumes to enable continuum, infinite-volume extrapolations [30, 34, 35].
Vector form factor results show encouraging consistency between LQCD and
experimental determinations [36, 29, 34]. Excited-state effects involving
$N\pi$ states have been identified as a significant source of systematic
uncertainty in LQCD calculations of axial form factors [37, 38, 28, 30, 34,
39]. More sophisticated analysis strategies employed in recent calculations
[28, 30, 34], as well as variational methods used to study $\pi\pi$ [40, 41,
42, 43, 44], $N\pi$ [45, 46, 47, 48, 49], and $NN$ [50, 51, 52] scattering and
$N\pi$ transition form factors [39] that enable $N\pi$ and other excited-state
effects to be explicitly subtracted from future nucleon elastic form factor
calculations, provide paths towards better quantifying and reducing these
challenging systematic uncertainties. Although significant future progress is
expected, LQCD calculations of nucleon axial form factors have progressed to a
point where it is important to understand the phenomenological impact of
current results and quantitatively establish what form factor precision is
required to achieve the cross-section precision needs of current and future
neutrino oscillation experiments [53, 54, 55].

Current LQCD results predict a significantly larger axial form factor at
$Q^{2}\sim 1$ GeV2 than determinations from deuterium bubble chamber data.111A
similar trend has recently been obtained within continuum Schwinger function
methods [56, 57]. LQCD form factor results were incorporated into the GENIE
[58] neutrino event generator and shown to lead to significant differences in
neutrino-nucleus cross sections in Ref. [54]. The interplay of these form
factor changes with other aspects of GENIE remains to be studied, and in
particular the “empirical MEC model” of two-body currents in GENIE and tuned
versions of this model used in neutrino oscillation experiments [59, 60, 61]
are obtained by fitting cross-section predictions to data while assuming a
particular axial form factor model for quasi-elastic contributions. It is
important to understand axial form factor effects on the event generators used
to analyze current experiments, but using such a data-driven approach to two-
body current contributions makes it difficult to disentangle the effects and
uncertainties associated with nucleon form factors from those associated with
two-body currents and other nuclear effects.

The theoretically well-defined separation between one- and two-body currents
in the GFMC and SF results of this work, in conjunction with the study of
nuclear model dependence enabled by comparing two realistic many-body methods,
enables robust quantification of the effects and uncertainties of nucleon
axial form factors and other single-nucleon inputs to neutrino-nucleus cross
sections. Parameterizing the nucleon axial form factor using the model-
independent $z$ expansion [62, 63, 64, 65] allows uncertainty quantification
to be performed without introducing any dependence on assumptions about the
shape of the axial form factor. Relations between neutrino-nucleus cross
section uncertainties and $z$ expansion parameter uncertainties are determined
below and used to construct quantitative precision targets for future LQCD
calculations of nucleon axial form factors. The effects of
$N\rightarrow\Delta$ transition form factors uncertainties are also studied.

This paper is organized as follows. Section II introduces the formalism
necessary to compute the inclusive neutrino-nucleus cross section. Section III
is dedicated to the determination of the axial form factor and its
uncertainties estimation. Results are presented in Section IV, while Section V
provides concluding remarks and discusses the outlook for future work.

##  II Methods

In the one-boson exchange approximation, the neutrino-nucleus double
differential cross sections can be written in the form

$\displaystyle\left(\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}\right)_{\nu/\bar{\nu}}$
$\displaystyle=\frac{G^{2}}{2\pi}\frac{k^{\prime}}{2E_{\nu}}[L_{CC}R_{CC}+2L_{CL}R_{CL}$
(1) $\displaystyle+L_{LL}R_{LL}+L_{T}R_{T}\pm 2L_{T^{\prime}}R_{T^{\prime}}],$

where the $\pm$ sign corresponds to $\nu(\bar{\nu})$ scattering. We take
$G=G_{F}\cos\theta_{c}$, with $G_{F}=1.1664\times 10^{-5}\,\rm{GeV^{-2}}$ [66]
and $\cos\theta_{c}=0.9740$ [67] for charged current processes studied here.
The $L_{i}$ factors only depend upon the kinematics of the initial and final
state leptons, while the electroweak response functions $R_{i}$ are defined as
linear combinations of different components of the hadronic response tensor
$R^{\mu\nu}$

$R^{\mu\nu}=\sum_{f}\langle 0|J^{\mu\dagger}|f\rangle\langle
f|J^{\nu}|0\rangle\delta(E_{0}+\omega-E_{f}).$ (2)

where $|0\rangle$ is the nuclear ground state, $|f\rangle$ are all possible
final states of the $A$-nucleon system and $J^{\mu}$ is the nuclear current
operator. The response tensor contains all the information on the structure of
the nuclear target, defined in terms of a sum over all transitions from the
ground state to any final state, including states with additional hadrons. In
this work we consider the one- and two-body contributions to the nuclear
current operator

$J^{\mu}=\sum_{i}j^{\mu}_{i}+\sum_{j>i}j^{\mu}_{ij}.$ (3)

and we will exclusively focus on final states involving only nucleons,
ignoring single-nucleon excitation processes. Exact expressions for $L_{i}$
and $R_{i}$ can be found in Ref. [68]. Computing the electroweak response
functions in the energy regime relevant to oscillation experiments is a
highly-nontrivial task. Existing approaches differ based on approximation
schemes and kinematic regimes of applicability. In this work, we will consider
two different methods that model the initial target state in a similar
fashion, but differ in the treatment of the interaction vertex and final-state
interactions: the Green’s function Monte Carlo, and the spectral function
approaches.

###  II.1 Green’s Function Monte Carlo

The response function in Eq. (2) is expressed in the energy domain. It can be
related to a matrix element of electroweak currents in the time domain,

$\displaystyle R_{\alpha\beta}(\mathbf{q},\omega)$ $\displaystyle=$
$\displaystyle\int dte^{i(\omega+E_{0})t}\ \langle 0|J_{\alpha}^{\dagger}\
e^{-iHt}\ J_{\beta}|0\rangle,$ (4)

where a completeness relation has been inserted to carry out the sum over the
possible final states of the $A$-nucleon system. The non-relativistic nuclear
Hamiltonian $H$ used in these calculations is comprised of the Argonne
$v_{18}$ [69] (AV18) nucleon-nucleon potentials plus the Illinois-7 [70] (IL7)
three-nucleon force. The Green’s function Monte Carlo method uses an
imaginary-time projection to distill information on the energy dependence of
the response functions. In particular, the imaginary-time response is obtained
by replacing the real time entering Eq. (4) with the imaginary time
$t\to-i\tau$. The Laplace transform of the energy-dependent response, called
the Euclidean response function, is defined by

$E_{\alpha\beta}(\mathbf{q},\tau)=\int_{\omega_{\rm
th}}^{\infty}d\omega\,e^{-\omega\tau}R_{\alpha\beta}(\mathbf{q},\omega)\,,$
(5)

where $\omega_{\rm th}$ is the inelastic threshold. Bayesian techniques, most
notably maximum entropy [71, 72], can be used to retrieve the energy
dependence of the response functions from their Laplace transform. Following
this strategy, GFMC calculations have been successfully performed to obtain
the inclusive electroweak response function of light nuclei while fully
retaining the complexity of many-body correlations and associated electroweak
currents [73, 74, 22]. Recently, algorithms based on artificial neural
networks have been developed to invert the Laplace transform [75] and show
better accuracy that maximum entropy, especially in the low-energy transfer
region.

The five response functions entering the CC cross section have been calculated
with GFMC methods for momentum transfers in the range (100–700) MeV in steps
of 100 MeV. In order to compute the flux-averaged cross sections presented in
Sec. IV that involve larger momentum transfers, we capitalize on the scaling
features of the GFMC response functions discussed in Refs. [76, 22, 77] to
interpolate and safely extrapolate them at momentum transfers
$|\mathbf{q}|>700$ MeV. The charge-changing weak current is the sum of vector
and axial components

$\displaystyle J^{\mu}$ $\displaystyle=J^{\mu}_{V}-J^{\mu}_{A},$
$\displaystyle J_{\pm}$ $\displaystyle=J_{x}\pm iJ_{y}\,.$ (6)

These operators comprise both one- and two-body contributions, whose
expressions are reported in Ref. [68]. Here, we focus on the axial and pseudo-
scalar one-body contributions given by

$\displaystyle j^{0}_{A}$
$\displaystyle=-\frac{F_{A}}{2m_{N}}\tau_{\pm}\bm{\sigma}\cdot\\{{\bf
k},e^{i{\bf q}\cdot{\bf r}}\\}$ $\displaystyle{\bf j}_{A}$
$\displaystyle=-F_{A}\tau_{\pm}\Big{[}{\bm{\sigma}}e^{i{\bf q}\cdot{\bf
r}}-\frac{1}{4m_{N}^{2}}\Big{(}{\bm{\sigma}}\\{{\bf k}^{2},e^{i{\bf
q}\cdot{\bf r}}\\}-\\{({\bm{\sigma}}\cdot{\bf k}){\bf k},e^{i{\bf q}\cdot{\bf
r}}\\}$ $\displaystyle-\frac{1}{2}{\bm{\sigma}\cdot{\bf q}}\\{{\bf k},e^{i{\bf
q}\cdot{\bf r}}\\}-\frac{1}{2}{\bf q}\\{({\bm{\sigma}}\cdot{\bf k}),e^{i{\bf
q}\cdot{\bf r}}\\}+i{\bf q}\times{\bf k}e^{i{\bf q}\cdot{\bf
r}}\Big{)}\Big{]}$ $\displaystyle j^{\mu}_{P}$
$\displaystyle=\frac{F_{P}}{2m_{N}^{2}}\tau_{\pm}q^{\mu}{\bm{\sigma}}\cdot{\bf
q}e^{i{\bf q}\cdot{\bf r}}\,,$ (7)

where $m_{N}$ is the nucleon mass and $\\{\cdot,\cdot\\}$ denotes the
anticommutator. The parametrization of the axial form factor $F_{A}$ has
historically involved models with and without over-constraining theoretical
assumptions and will be a focus of this work in Sec. III.1. For the pseudo-
scalar form factor, PCAC and pion-pole dominance relations valid in leading
order chiral perturbation theory [78, 79] can be used to relate $F_{P}$ to
$F_{A}$ as

$F_{P}=\frac{2F_{A}m_{N}^{2}}{m^{2}_{\pi}+Q^{2}},$ (8)

where $Q^{2}=-q^{2}=-(q^{0})^{2}+\mathbf{q}^{2}$ and $m_{\pi}$ is the pion
mass. This relation is consistent with the pseudoscalar form factor
constraints extracted from precise measurements of the muon-capture rate on
hydrogen and 3He [80]. It is also consistent with current LQCD results [30,
34]. Precise future LQCD calculations will further test Eq. (8) and provide
independent determinations of $F_{P}$ without theoretical assumptions on its
relation to $F_{A}$.

The isovector two-body currents can be separated into model-dependent and
-independent terms. Those associated with pion exchange are denoted as model
independent; they are constrained by the continuity equation and do not
contain any free parameters, since they are determined directly from the
nucleon-nucleon interaction. The expressions adopted in this work and reported
in Refs. [68, 81] are consistent with the semi-phenomenological AV18
interaction. The diagrams associated with the intermediate excitation of a
$\Delta$ isobar are purely transverse and generally referred to as model
dependent. Their expression can be found in Refs. [68, 81], it is important to
mention that in order to be used in the GFMC calculations, the static $\Delta$
approximation has to be adopted, i.e. the kinetic-energy contributions in the
denominator of the $\Delta$ propagator are neglected. The impact of this
approximations will be further discussed when comparing the GFMC and SF
results. Among the axial two-body current operators, the leading terms of
pionic range are those associated with excitation of $\Delta$-isobar
resonances, despite being less relevant the two-body operators associated with
axial $\pi NN$ contact interactions are also accounted for, we refer the
reader to Refs [68, 81] for a more detailed discussion of these contributions.

###  II.2 Extended Factorization Scheme

At large values of momentum transfer ($|\textbf{q}|\gtrsim 400$ MeV), a
factorization scheme can be used in which the neutrino-nucleus scattering is
approximated as an incoherent sum of scatterings with individual nucleons, and
the struck nucleon system is decoupled from the rest of the final state
spectator system. In the quasielastic region, the dominant reaction mechanism
is single nucleon knockout. In this case, the final state is factorized
according to

$|f\rangle=|\textbf{p}^{\prime}\rangle\otimes|\Psi^{A-1}_{f},\textbf{p}_{A-1}\rangle,$
(9)

where $|\textbf{p}^{\prime}\rangle$ is the final state nucleon produced at the
vertex, assumed to be in a plane wave state, and
$|\Psi^{A-1}_{f},\textbf{p}_{A-1}\rangle$ describes the residual system,
carrying momentum $\textbf{p}_{A-1}$. Inserting this factorization ansatz as
well as a single-nucleon completeness relation gives the matrix element of the
1 body current operator as

$\langle
f|j^{\mu}|0\rangle\rightarrow\sum_{k}[\langle\Psi^{A-1}_{f}|\otimes\langle
k|]|0\rangle\langle p|\sum_{i}j^{\mu}_{i}|k\rangle,$ (10)

where $\textbf{p}=\textbf{q}+\textbf{k}$. This first piece of the matrix
element explicitly does not depend on the momentum transfer and so can be
computed using techniques in nuclear many body theory. The second piece can be
straightforwardly computed once the currents $j^{\mu}_{i}$ are specified as
the single nucleon states are just free Dirac spinors. Substituting the last
equation into Eq. (2), and exploiting momentum conservation at the single
nucleon vertex, allows us to rewrite the one body contribution to the response
tensor as

$\displaystyle R_{1b}^{\mu\nu}(\textbf{q},\omega)=$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{3}}dEP_{h}(\textbf{k},E)\frac{m^{2}_{N}}{e(\textbf{k})e(\textbf{k}+\textbf{q})}$
(11) $\displaystyle\times\sum_{i}\langle
k|j_{i}^{\mu\dagger}|k+q\rangle\langle k+q|j_{i}^{\nu}|k\rangle$
$\displaystyle\times\delta(\tilde{\omega}+e(\textbf{k})-e(\textbf{p})),$

where $e(\mathbf{k})=\sqrt{m_{N}^{2}+\mathbf{k}^{2}}$. The factors
$m_{N}/e(\textbf{k})$ and $m_{N}/e(\textbf{k}+\textbf{q})$ are included to
account for the covariant normalization of the four spinors in the matrix
elements of the relativistic current. The energy transfer has been replaced by
$\tilde{\omega}=\omega-m_{N}+E-e(\textbf{k})$ to account for some of the
initial energy transfer going into the residual nuclear system. Finally, the
calculation of the one-nucleon spectral function $P_{h}(\textbf{k},E)$
provides the probability of removing a nucleon with momentum k and leaving the
residual nucleus with an excitation energy $E$. Its derivation using QMC
techniques is discussed in Sec. II.3. The relativistic current used in Eq.
(11) is given as a sum of vector and axial terms which can be written as

$\displaystyle
j^{\mu}_{V}=F_{1}\gamma^{\mu}+i\sigma^{\mu\nu}q_{\nu}\frac{F_{2}}{2m_{N}}$
(12) $\displaystyle
j^{\mu}_{A}=\gamma^{\mu}\gamma^{5}F_{A}+q^{\mu}\gamma^{5}\frac{F_{P}}{m_{N}}.$

In the above, $F_{1}$ and $F_{2}$ are usual isovector and isoscalar form
factors which themselves are functions of the proton and neutron electric and
magnetic form factors. These are highly constrained by electron scattering and
we adopt the Kelly parameterization in this work [82]. The expression of the
pseudo scalar form factor in terms of the axial one is the same of Eq. (8).
The details on the different parameterizations adopted for $F_{A}$ are given
in Secs. III. To treat amplitudes involving two-nucleon currents, the
factorization ansatz of Eq. (9) can be generalized as

$\displaystyle|\psi_{f}^{A}\rangle\rightarrow|pp^{\prime}\rangle_{a}\otimes|\psi_{f}^{A-2}\rangle\,.$
(13)

where $|p\,p^{\prime}\rangle_{a}=|p\,p^{\prime}\rangle-|p^{\prime}\,p\rangle$
is the anti-symmetrized state of two-plane waves with momentum $p$ and
$p^{\prime}$. Following the work presented in Refs. [83, 84, 23], the pure
two-body current component of the response tensor can be written as

$\displaystyle R^{\mu\nu}_{\rm 2b}({\bf q},\omega)=\frac{V}{2}\int
dE\frac{d^{3}k}{(2\pi)^{3}}\frac{d^{3}k^{\prime}}{(2\pi)^{3}}\frac{d^{3}p}{(2\pi)^{3}}\frac{m_{N}^{4}}{e({\bf
k})e({\bf k^{\prime}})e({\bf p})e({\bf p^{\prime}})}$
$\displaystyle\qquad\times P_{h}^{\rm NM}({\bf k},{\bf
k}^{\prime},E)\sum_{ij}\,\langle
k\,k^{\prime}|{j_{ij}^{\mu}}^{\dagger}|p\,p^{\prime}\rangle_{a}\langle
p\,p^{\prime}|j_{ij}^{\nu}|k\,k^{\prime}\rangle$
$\displaystyle\qquad\times\delta(\omega-E+2m_{N}-e(\mathbf{p})-e(\mathbf{p}^{\prime}))\,.$
(14)

In the above equation, the normalization volume for the nuclear wave functions
$V=\rho/A$ with $\rho=3\pi^{2}k_{F}^{3}/2$ depends on the Fermi momentum of
the nucleus, which for 12C is taken to be $k_{F}=225$ MeV. The factor $1/2$
accounts for the double counting that occurs in this notation (the product of
the two direct terms is equal to the one of the two exchange terms).
Differently from Refs. [23], the two-nucleon spectral function adopted in this
work is not factorized into the product of two one-nucleon spectral functions,
see Sec. II.3 for a detailed discussion.

The two-body CC operator is given by the sum of four distinct interaction
mechanisms, namely the pion in flight, seagull, pion-pole, and $\Delta$
excitations [85, 86]

$\displaystyle j^{\mu}_{\rm CC}=(j^{\mu}_{\rm pif})_{\rm CC}+(j^{\mu}_{\rm
sea})_{\rm CC}+(j^{\mu}_{\rm pole})_{\rm CC}+(j^{\mu}_{\Delta})_{\rm CC}.$
(15)

Figure 1: Feynman diagrams describing of the first two contributions to the
two-body currents associated with $\Delta$-excitation processes. Solid, thick
green, and dashed lines correspond to nucleons, deltas, pions, respectively.
The wavy line represents the vector boson.

Detailed expressions for the first four terms of Eq. (15) can be found in
Refs. [86, 23]. Below, we only report the two-body current terms involving a
$\Delta$-resonance in the intermediate state, as we find them to be the
dominant contribution. Because of the purely transverse nature of this
current, the form of its vector component is not subject to current-
conservation constraints and its expression is largely model dependent, as
discussed in the previous section. The current operator can be written as [85,
86]:

$\displaystyle(j^{\mu}_{\Delta})_{\rm CC}$
$\displaystyle=\frac{3}{2}\frac{f_{\pi
NN}f^{\ast}}{m^{2}_{\pi}}\Bigg{\\{}\Bigg{[}\Big{(}-\frac{2}{3}\tau^{(2)}_{\pm}+\frac{(\tau^{(1)}\times\tau^{(2)})_{\pm}}{3}\Big{)}$
$\displaystyle\times F_{\pi NN}(k^{\prime}_{\pi})F_{\pi
N\Delta}(k^{\prime}_{\pi})(j^{\mu}_{a})_{(1)}$
$\displaystyle-\Big{(}\frac{2}{3}\tau^{(2)}_{\pm}+\frac{(\tau^{(1)}\times\tau^{(2)})_{\pm}}{3}\Big{)}_{\pm}\Big{)}$
$\displaystyle\times F_{\pi NN}(k^{\prime}_{\pi})F_{\pi
N\Delta}(k^{\prime}_{\pi})(j^{\mu}_{b})_{(1)}\Bigg{]}\Pi(k^{\prime}_{\pi})_{(2)}+(1\leftrightarrow
2)\Bigg{\\}}$ (16)

where $k^{\prime}$ and $p^{\prime}$ are the initial and final momentum of the
second nucleon, respectively, while $k^{\prime}_{\pi}=p^{\prime}-k^{\prime}$
is the momentum of the $\pi$ exchanged in the two depicted diagrams of Fig. 1,
$f^{\ast}$=2.14, and

$\displaystyle\Pi(k_{\pi})=\frac{\gamma_{5}\not{k}_{\pi}}{k_{\pi}^{2}-m^{2}_{\pi}}\
,$ (17) $\displaystyle F_{\pi N\Delta}(k_{\pi})=\frac{\Lambda^{2}_{\pi
N\Delta}}{\Lambda^{2}_{\pi N\Delta}-k_{\pi}^{2}}\ ,$ (18) $\displaystyle
F_{\pi
NN}(k_{\pi})=\frac{\Lambda_{\pi}^{2}-m_{\pi}^{2}}{\Lambda_{\pi}^{2}-k_{\pi}^{2}}\
,$ (19)

with $\Lambda_{\pi N\Delta}=1150$ MeV and $\Lambda_{\pi}=1300$ MeV. The term
$\tau_{\pm}=(\tau_{x}\pm i\tau_{y})/2$ is the isospin raising/lowering
operator. In Eq. (16), $j^{\mu}_{a}$ and $j^{\mu}_{b}$ denote the
$N\rightarrow\Delta$ transition vertices of diagram (a) and (b) of Fig. 1,
respectively. The expression of $j^{\mu}_{a}$ is given by

$\displaystyle j^{\mu}_{a}$
$\displaystyle=(j^{\mu}_{a})_{V}+(j^{\mu}_{a})_{A}\ ,$
$\displaystyle(j^{\mu}_{a})_{V}$
$\displaystyle=(k_{\pi}^{\prime})^{\alpha}G_{\alpha\beta}(p_{\Delta})\Big{[}\frac{C_{3}^{V}}{m_{N}}\Big{(}g^{\beta\mu}\not{q}-q^{\beta}\gamma^{\mu}\Big{)}$
$\displaystyle+\frac{C_{4}^{V}}{m_{N}^{2}}\Big{(}g^{\beta\mu}q\cdot
p_{\Delta}-q^{\beta}p_{\Delta}^{\mu}\Big{)}$
$\displaystyle+\frac{C_{5}^{V}}{m_{N}^{2}}\Big{(}g^{\beta\mu}q\cdot
k-q^{\beta}k^{\mu}+C_{6}^{V}g^{\beta\mu}\Big{)}\Big{]}\gamma_{5}\ ,$
$\displaystyle(j^{\mu}_{a})_{A}$
$\displaystyle=(k_{\pi}^{\prime})^{\alpha}G_{\alpha\beta}(p_{\Delta})\Big{[}\frac{C_{3}^{A}}{m_{N}}\Big{(}g^{\beta\mu}\not{q}-q^{\beta}\gamma^{\mu}\Big{)}$
$\displaystyle+\frac{C_{4}^{A}}{m_{N}^{2}}\Big{(}g^{\beta\mu}q\cdot
p_{\Delta}-q^{\beta}p_{\Delta}^{\mu}\Big{)}$
$\displaystyle+C_{5}^{A}g^{\beta\mu}+\frac{C_{6}^{A}}{m_{N}^{2}}q^{\mu}q^{\alpha}\Big{]},$
(20)

where $k$ is the momentum of the initial nucleon which absorbs the incoming
momentum $\tilde{q}$ and $p_{\Delta}=\tilde{q}+k$, yielding
$p^{0}_{\Delta}=e({\bf k})+\tilde{\omega}$. We introduced
$\tilde{q}=(\tilde{\omega},{\bf q})$ to account for the fact that the initial
nucleons are off-shell. A similar definition can be written down for
$j^{\mu}_{b}$; more details are reported in Ref. [24, 23]. For $C_{3}^{V}$ we
adopted the model of Ref. [87], yielding

$\displaystyle
C_{3}^{V}=\frac{2.13}{(1-q^{2}/M_{V}^{2})^{2}}\frac{1}{1-q^{2}/(4M_{V}^{2})},$
(21)

with $M_{V}=0.84$ GeV. Following the discussion of Ref. [86], we neglected the
terms $C_{4}^{V}$ and $C_{5}^{V}$ which are expected to be suppressed by
$\mathcal{O}(k/m_{N})$, while $C_{6}^{V}=0$ by conservation of the vector
current. However, it is worth mentioning that including these terms in the
current operator would not pose any conceptual difficulty. To be consistent,
in the axial part we only retain the leading contribution of Eq. (20), which
is the term proportional to $C_{5}^{A}$ defined as [88]

$\displaystyle
C_{5}^{A}=\frac{1.2}{(1-q^{2}/M_{A\Delta})^{2}}\times\frac{1}{1-q^{2}/(3M_{A\Delta})^{2})},$
(22)

with $M_{A\Delta}=1.05$ GeV.

The Rarita-Schwinger propagator

$\displaystyle
G^{\alpha\beta}(p_{\Delta})=\frac{P^{\alpha\beta}(p_{\Delta})}{p^{2}_{\Delta}-M_{\Delta}^{2}},$
(23)

is proportional to the spin 3/2 projection operator
$P^{\alpha\beta}(p_{\Delta})$. In order to account for the possible decay of
the $\Delta$ into a physical $\pi N$, we replace $M_{\Delta}\rightarrow
M_{\Delta}-i\Gamma(p_{\Delta})/2$ [89, 90] where the last term is the energy
dependent decay width given by

$\displaystyle\Gamma(p_{\Delta})$ $\displaystyle=\frac{(4f_{\pi
N\Delta})^{2}}{12\pi
m_{\pi}^{2}}\frac{|\mathbf{d}|^{3}}{\sqrt{s}}(m_{N}+E_{d})R(\mathbf{r}^{2})\,.$
(24)

In the above equation, $(4f_{\pi N\Delta})^{2}/(4\pi)=0.38$,
$s=p_{\Delta}^{2}$ is the invariant mass, $\mathbf{d}$ is the decay three-
momentum in the $\pi N$ center of mass frame, such that

$|\mathbf{d}|^{2}=\frac{1}{4s}[s-(m_{N}+m_{\pi})^{2}][s-(m_{N}-m_{\pi})^{2}]\,$
(25)

and $E_{d}=\sqrt{m_{N}^{2}+\mathbf{d}^{2}}$ is the associated energy. The
additional factor

$R(\mathbf{r}^{2})=\left(\frac{\Lambda_{R}^{2}}{\Lambda_{R}^{2}-\mathbf{r}^{2}}\right),$
(26)

depending on the $\pi N$ three-momentum $\mathbf{r}$, with
$\mathbf{r}^{2}=(E_{d}-\sqrt{m_{\pi}^{2}+\mathbf{d}^{2}})^{2}-4\mathbf{d}^{2}$
and $\Lambda_{R}^{2}=0.95\,m_{N}^{2}$, is introduced to improve the
description of the experimental phase-shift $\delta_{33}$ [89]. The medium
effects on the $\Delta$ propagator are accounted for by modifying the decay
width as

$\displaystyle\Gamma_{\Delta}(p_{\Delta})\to\Gamma_{\Delta}(p_{\Delta})-2\rm{Im}[U_{\Delta}(p_{\Delta},\rho=\rho_{0})],$
(27)

where $U_{\Delta}$ is a density dependent potential obtained from a Bruckner-
Hartree-Fock calculation using a coupled-channel $NN\oplus N\Delta\oplus\pi
NN$ model [91, 92, 93, 94] and we fixed the density at the nuclear saturation
value $\rho_{0}$ =0.16 fm3. For a detailed analysis of medium effects in the
MEC contribution for electron-nucleus scattering see Ref. [24].

###  II.3 Spectral Function

The spectral function of a nucleon with isospin $\tau_{k}=p,n$ and momentum
${\bf k}$ can be written as

$\displaystyle P_{\tau_{k}}(\mathbf{k},E)$ $\displaystyle=\sum_{n}|\langle
0|[|k\rangle\,|\Psi_{n}^{A-1}\rangle]|^{2}\delta(E+E_{0}-E_{n}^{A-1})\,.$ (28)

where $|k\rangle$ is the single-nucleon state, $|0\rangle$ is the ground state
of the Hamiltonian with energy $E_{0}$, while $|\Psi_{n}^{A-1}\rangle$ and
$E_{n}^{A-1}$ are the energy eigenstates and eigenvalues of the remnant
nucleus with $(A-1)$ particles. The momentum distribution of the initial
nucleon is obtained by integrating the spectral function over the removal
energy

$n_{\tau_{k}}(k)=\int dEP_{\tau_{k}}(\mathbf{k},E)\,.$ (29)

We can rewrite the spectral function as a sum of a mean field (MF) and a
correlation (corr) term. For simplicity, let us focus on the $\tau_{k}=p$
case; the MF contribution of the proton spectral function is obtained by
considering only bound $A-1$ states of the remnant nucleus

$\displaystyle P_{p}^{\rm MF}(\mathbf{k},E)$ $\displaystyle=\sum_{n}|\langle
0^{A}|[|k\rangle\otimes|\Psi_{n}^{A-1}\rangle]|^{2}$
$\displaystyle\times\delta\Big{(}E-B^{A}_{0}+B_{n}^{A-1}-\frac{k^{2}}{2m^{A-1}_{n}}\Big{)}\,,$
(30)

where $B^{A}_{0}$ and $B_{n}^{A-1}$ are the binding energies of the initial
and the bound $A-1$ spectator nucleus with mass and $m_{n}^{A-1}$. The
momentum-space overlaps $\langle
0^{A}|[|k\rangle\otimes|\Psi_{n}^{A-1}\rangle$ pertaining to the p-shell
contributions are computed by Fourier transforming the variational Monte Carlo
(VMC) radial overlaps for the transitions [95]:

${}^{12}{\rm C}(0^{+})$ $\displaystyle\rightarrow^{11}{\rm B}(3/2^{-})+p$
${}^{12}{\rm C}(0^{+})$ $\displaystyle\rightarrow^{11}{\rm B}(1/2^{-})+p$
${}^{12}{\rm C}(0^{+})$ $\displaystyle\rightarrow^{11}{\rm
B}(3/2^{-})^{\ast}+p\,.$

The calculation of the s-shell mean-field contribution involves non trivial
difficulties for the VMC method, as it would require to evaluate the
spectroscopic overlaps for the transitions to all the possible excited states
of ${}^{11}{\rm B}$ with $J^{P}=(1/2^{+})$. To overcome this limitation, we
used the VMC overlap associated with the ${}^{4}{\rm
He}(0^{+})\rightarrow{}^{3}{\rm H}(1/2^{+})+p$ transition and applied minimal
changes to the quenching factor which is needed to reproduce the integral of
the momentum distribution up to $k_{F}=1.15$ fm-1. More details about the
adopted procedure are discussed in Ref. [26].

The correlation contribution to the SF is given by

$\displaystyle P_{p}^{\rm corr}$
$\displaystyle(\mathbf{k},E)=\sum_{n}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}|\langle
0^{A}|[|k\rangle\,|k^{\prime}\rangle\,|\Psi_{n}^{A-2}\rangle]|^{2}$
$\displaystyle\times\delta(E+E_{0}-e(\mathbf{k}^{\prime})-E_{n}^{A-2})\,$
$\displaystyle=\mathcal{N}_{p}\sum_{\tau_{k^{\prime}}=p,n}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}}\Big{[}n_{p,\tau_{k^{\prime}}}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle\qquad\times\delta\Big{(}E-B_{0}-e(\mathbf{k}^{\prime})+\bar{B}_{A-2}-\frac{(\mathbf{k}+\mathbf{k}^{\prime})^{2}}{2m_{A-2}}\Big{)}\Big{]}\,,$
(31)

to derive the last expression, we used a completeness relation and assumed
that the $(A-2)$-nucleon binding energy is narrowly distributed around a
central value $\bar{B}_{A-2}$. The mass of the recoiling $A-2$ system is
denoted by $m_{A-2}$ and $\mathcal{N}_{p}$ is an appropriate normalization
factor. We started from the VMC two-nucleon momentum distribution
$n_{\tau_{k},\tau_{k^{\prime}}}(\mathbf{k},\mathbf{k}^{\prime})$ of Ref. [96],
but in order to isolate the contribution of short-range correlated nucleons we
performed cuts in the relative momentum of the pairs, requiring that the
overall normalization and shape of the one-nucleon momentum distributions are
correctly recovered.

In this work, we only consider the MF contribution to the two-nucleon spectral
function, e.g. we neglect contributions where more than two nucleons are
emitted. We write $P_{\tau_{k},\tau_{k}^{\prime}}^{\rm MF}$ as

$\displaystyle P_{\tau_{k},\tau_{k}^{\prime}}^{\rm
MF}(\mathbf{k},\mathbf{k}^{\prime},E)$
$\displaystyle=n_{\tau_{k},\tau_{k^{\prime}}}(\mathbf{k},\mathbf{k}^{\prime})$
$\displaystyle\times\delta\Big{(}E-B_{0}+\bar{B}_{A-2}-\frac{\mathbf{K}^{2}}{2m_{A-2}}\Big{)}\,.$
(32)

where $\mathbf{K}=\mathbf{k}+\mathbf{k}^{\prime}$ is the total momentum of the
pair. It has been argued that the strong isospin dependence of the two-nucleon
momentum distribution, supported by experimental data, persists for nuclei
even larger than 12C [97, 98, 99, 100]. In Refs. [23, 24]
$P(\mathbf{k},\mathbf{k}^{\prime},E)$ has been written as the product of two
one-nucleon spectral functions assuming that long-range correlations are
absent. In the present work, we go beyond this approximation and adopt a novel
definition of the two-nucleon spectral function that retains correlations
between the two struck particles. A detailed comparison of the results
obtained with these different parametrizations of the two-nucleon spectral
functions will be carried out in a future work and validated against electron
scattering data.

##  III Nucleon axial form factor effects and uncertainties

Imprecise knowledge of the nucleon axial form factor, $F_{A}$, might lead to
significant systematic uncertainties in neutrino-nucleus cross-section
calculations. This is in contrast with the nucleon vector form factors in Eq.
(12), which can be precisely determined using high-statistics electron
scattering experiments [101, 102, 103, 104, 105]. Current experimental
constraints on nucleon axial form factors come from beta decay measurements,
neutrino scattering on nuclear targets, and pion electroproduction [4, 5, 6,
7, 8, 10, 11, 12, 64, 65]. These give weak constraints on $F_{A}$ in
comparison to the vector form factors, as beta decay is only sensitive the
absolute normalization $g_{A}=F_{A}(0)$ and neutrino scattering and pion
electroproduction experiments are limited by both statistics and nuclear
modeling uncertainties.

As LQCD calculations of nucleon form factors mature [27, 28, 30, 31, 32, 33,
34, 35], it becomes increasingly important to quantify the level of axial form
factor precision required to achieve a given level of neutrino-nucleus cross-
section accuracy. This is challenging because axial form factor effects on
flux-averaged neutrino-nucleus cross sections can be difficult to disentangle
from nuclear effects such as two-body currents, as evident for example in the
differences between theoretical descriptions of MiniBooNE data with either an
unexpectedly slow falloff of the axial form factor with increasing momentum
transfer [13] or with larger than anticipated contributions from two-body
current effects [15, 16, 17, 18, 19, 20, 21]. This ambiguity between one- and
two-body current effects on flux-averaged cross sections makes it essential to
quantify the role of the nucleon axial form factor in neutrino-nucleus cross-
section calculations using nuclear effective theories that provide a
consistent theoretical decomposition between one- and two-body current
contributions. The remainder of this section discusses how to quantify nucleon
axial form factor effects on neutrino-nucleus cross-section calculations based
on the model-independent $z$ expansion and how to estimate nucleon axial form
faction precision needs using the GFMC and spectral function methods discussed
above.

###  III.1 Parametrization

Historically a dipole parametrization has often been used for the axial form
factor

$F_{A}(Q^{2})=\frac{g_{A}}{(1+Q^{2}/M_{A}^{2})^{2}},$ (33)

where $g_{A}=1.2723(23)$ has been measured from neutron beta decay [106], and
$M_{A}=1.014\pm 0.014$ GeV [107]. However, this one-parameter form is not
expressive enough to describe the shape of the axial form factor predicted by
QCD. It has been demonstrated in Refs. [64, 65] that assuming that the dipole
parameterization is valid when fitting to experimental results can lead to
form factor fits with uncertainties that are underestimated by a factor of
$\sim 5$ in comparison to those determined using fits based off a model-
independent $z$ expansion. Axial form factors in QCD are analytic functions of
$Q^{2}$ for $Q^{2}=-t>-t_{c}$, where $t_{c}$ is the location of the
$t$-channel cut, which enables an analytic function $z(Q^{2})$ to be defined
as [62, 63, 64]

$z(Q^{2})=\frac{\sqrt{t_{c}+Q^{2}}-\sqrt{t_{c}-t_{0}}}{\sqrt{t_{c}+Q^{2}}+\sqrt{t_{c}-t_{0}}},$
(34)

where $t_{0}$ is an arbitrary parameter whose choice is discussed in Sec.
III.2 below. For $F_{A}$ the cut starts at $t_{c}=9m_{\pi}^{2}$. Because
$|z|<1$, the axial form factor can be expanded as a power series in $z(Q^{2})$
for the $Q^{2}>0$ domain of interest for neutrino-nucleus scattering,

$F_{A}(Q^{2})=\sum_{k=0}^{\infty}a_{k}\,z(Q^{2})^{k}\approx\sum_{k=0}^{k_{\rm
max}}a_{k}\,z(Q^{2})^{k},$ (35)

where the $z$ expansion coefficients $a_{k}$ include nucleon structure
information and $k_{\rm max}$ is a truncation parameter required to make the
number of expansion parameters finite. The parameter $a_{0}$ can be fixed by
the sum rule

$\sum_{k=0}^{k_{\rm max}}a_{k}z(0)^{k}=g_{A}.$ (36)

Constraints on the $a_{k}$ are also obtained by enforcing the correct large
$Q^{2}$ behavior of the axial form factor, which is predicted by perturbative
QCD to be $Q^{-4}$ up to logarithmic corrections [108]. This asymptotic
$Q^{-4}$ behavior can be enforced by demanding that $F_{A}(Q^{2})$ and its
first three derivatives with respect to $1/Q$ vanish for asymptotically large
$Q^{2}$, corresponding to $z=1$, which is equivalent to

$\left.\frac{d^{n}}{dz^{n}}F_{A}\right|_{z=1}=0\;\;\;;\;\;\;n=0,1,2,3,$ (37)

and therefore leads to the sum rules [109]

$\sum_{k=n}^{k_{\rm max}}\frac{k!}{(k-n)!}\ a_{k}=0\;\;\;;\;\;\;n=0,1,2,3.$
(38)

In practice, these constraints can be satisfied by first determining the
$k_{\rm max}$ and $a_{1},\ldots,a_{k_{\rm max}}$ preferred by a fit to data
(with $a_{0}$ either treated as an additional independent parameters or as
being fixed by the constraint Eq. (36)), and then replacing $k_{\rm max}$ with
$k_{\rm max}+4$ and solving for the four unconstrained coefficients using Eq.
(38). The remaining $a_{k}$ must then be fixed by information on the
$Q^{2}$-dependence of the axial form factor determined theoretically using
LQCD calculations or experimentally by fitting neutrino-nucleus scattering
and/or pion electroproduction data. The $z$ expansion can be used to provide a
model-independent definition of the dependence of neutrino-nucleus cross-
section uncertainties on nucleon axial form factor uncertainties. Any function
$\sigma(F_{A},X)$ that depends on the axial form factor, as well as any number
of additional independent form factors and parameters schematically denoted by
$X$, can be interpreted as a function of the axial form factor $z$ expansion
coefficients, $\sigma(a_{k},X)$. Neglecting correlations among the $a_{k}$ and
between $X$ and the $a_{k}$, the relative uncertainties of $\sigma$, $a_{k}$,
and $X$ can be expressed as

$\begin{split}\frac{\delta\sigma}{|\sigma|}&=\left|\frac{X}{\sigma}\right|\frac{\delta
X}{|X|}+\sum_{k}\left|\frac{\partial\sigma}{\partial
a_{k}}\right|\left|\frac{a_{k}}{\sigma}\right|\frac{\delta a_{k}}{|a_{k}|}\\\
&\geq\left|\frac{\partial\sigma}{\partial
a_{k}}\right|\left|\frac{a_{k}}{\sigma}\right|\frac{\delta
a_{k}}{|a_{k}|},\end{split}$ (39)

where the inequality holds for each $a_{k}$ individually by positivity of each
term in the sum. In particular, Eq. (39) can be applied to neutrino-nucleus
differential cross sections $\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}$ at the
values of $T_{\mu}$ that maximize the cross-section for a given
$\cos\theta_{\mu}$, denoted
$d\sigma_{P}\equiv\left.\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}\right|_{\rm
peak}$, to give

$\begin{split}\frac{\delta\left(d\sigma_{P}\right)}{d\sigma_{P}}&\geq
D_{k}(\theta_{\mu})\frac{\delta a_{k}}{|a_{k}|},\end{split}$ (40)

where

$D_{k}(\theta_{\mu})=\left|\frac{\partial d\sigma_{P}}{\partial
a_{k}}\right|\frac{\left|a_{k}\right|}{d\sigma_{P}}.$ (41)

A simple way to approximate the derivative appearing in $D_{k}(\theta)$ is to
compute the finite difference,

$\begin{split}\frac{\partial d\sigma_{P}}{\partial
a_{k}}\approx&\frac{1}{\delta a_{k}}\left[d\sigma_{P}(a_{k}+\delta
a_{k})-d\sigma_{P}(a_{k})\right],\end{split}$ (42)

where $\delta a_{k}$ is a small parameter and $d\sigma_{P}(a_{k}+\delta
a_{k})$ denotes the peak differential cross section calculated with $a_{k}$
replaced by $a_{k}+\delta a_{k}$ and $a_{i}$ held fixed for $i\neq k$. The
accuracy of the finite-difference approximation can be validated by performing
calculations with different values of $\delta a_{k}$ and verifying the
stability of the resulting derivative estimate.

Since neutrino mixing angle determinations at oscillation experiments are
sensitive to the relative heights of the cross-section peaks at the near and
far detectors, the relative uncertainty in
$\left.\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}\right|_{\rm peak}$ provides a
reasonable measure of the cross-section uncertainty entering neutrino
oscillation analyses.222Further information about the uncertainty needs of
neutrino oscillation experiments could be obtained by applying Eq. (39)
directly to the neutrino mass and mixing parameters derived from an
oscillation analysis with for example $\frac{\partial\delta_{CP}}{\partial
a_{k}}$ estimated in analogy to Eq. (42). This requires a detailed analysis
including for example neutrino energy reconstruction and detector efficiency
effects and is left to future work.  The relative uncertainty on $a_{k}$
required for axial form factor uncertainties to not saturate a given cross-
section uncertainty budget according to this measure is then given by the
desired cross-section uncertainty times $1/D_{k}(\theta)$.

###  III.2 Form Factor Results

Neutrino scattering experiments typically involve nuclear targets, and the
challenges of determining nucleon axial form factors from data arise both from
scarcity of precise experimental data and from the need to understand nuclear
effects that are present in data. Relatively theoretically clean
determinations can be made using neutrino-deuteron scattering data obtained
from bubble-chamber experiments at ANL, BNL, and FNAL in the 70s and 80s [4,
5, 6, 7, 8], which were recently analyzed in Ref. [65] using the $z$ expansion
to parameterize the nucleon axial form factor. It was found that uncertainties
on the axial form factor shape were significantly underestimated by previous
analysis that assumed an over-constraining dipole form factor shape, and in
particular the uncertainties on the axial radius are estimated to be five
times larger using fits to the $z$ expansion than using dipole fits. The
relation between the free nucleon and deuteron cross sections need for this
analysis is parameterized as a $Q^{2}$-dependent but $E_{\nu}$-independent
function taken from nuclear theory calculations [110, 68]. A 10% systematic
uncertainty is included in the analysis of Ref. [65] to account for nuclear
model uncertainties and other systematic corrections.

Figure 2: The nucleon axial form factor $F_{A}(Q^{2})$ determined using fits
to neutrino-deuteron scattering data using the model-independent $z$ expansion
from Ref. [65] (D2 Meyer et al.) are shown as a blue band in the top panel.
LQCD results are shown for comparison from Ref. [30] (LQCD Bali et al.,
green), Ref. [34] (LQCD Park et al., red) and Ref. [35] (LQCD Djukanovic et
al., purple). Bands show combined statistical and systematic uncertainties in
all cases, see the main text for more details. A dipole parameterization with
$M_{A}=1.0$ GeV and a $1.4\%$ uncertainty [107] is also shown for comparison
(black). The lower panel shows the absolute value of the difference between D2
Meyer et al. and LQCD Bali et al. results divided by their uncertainties added
in quadrature, denoted $\delta F_{A}/\sigma$; very similar results are
obtained using the other LQCD results.

Alternative strategies for constraining nucleon axial form factors without
nuclear uncertainties are to perform neutrino scattering experiments with
hydrogen targets, see Ref. [111] for a discussion of current proposals and
challenges, and to calculate nucleon axial form factors theoretically using
LQCD, as reviewed in Refs. [53, 54, 55]. There is a long history of LQCD
studies of nucleon form factors [112, 113] including recent and ongoing work
by multiple groups [27, 28, 29, 30, 31, 32, 33, 34, 35]. Challenging
systematic uncertainties arising from $N\pi$ excited-state contamination have
been identified in axial form factor calculations [37, 38, 28, 30, 34, 39],
and recent calculations explicitly accounting for these effects have shown
consistency with identities following from axial ward identities and the
assumption of negligible excited-state contamination [28, 30, 34]. Although
further study of these excited-state effects is needed [39], calculations have
now been performed using approximately physical quark masses and continuum,
infinite-volume extrapolations [30, 34, 35]. In each of these LQCD
calculations, axial form factors are parameterized using the $z$ expansion and
fit to nucleon axial-current matrix elements at several discrete values of
$Q^{2}$, though different strategies for parameterizing discretization effects
and other sources of systematic uncertainty are employed.333The uncertainties
reported in Ref. [34] are inflated by a factor of 3 to account for systematic
uncertainties arising from lattice spacing and quark mass effects [111].

The LQCD nucleon axial form factor results of Refs. [30, 34, 35] are shown in
Fig. 2 along with the form factor results determined from experimental
neutrino-deuteron scattering data in Ref. [65]. Fits were performed using
results with $Q^{2}\leq 1$ GeV2 in Refs. [65, 30, 34] and with $Q^{2}\leq 0.7$
GeV2 in Ref. [35] with the parameterization provided by the $z$ expansion used
to extrapolate form factor results to larger $Q^{2}$. Clear agreement between
different LQCD calculations can be seen. However, the LQCD axial form factor
results are 2-3$\sigma$ larger than the results of Ref. [65] for $Q^{2}\gtrsim
0.3\text{ GeV}^{2}$. The effects of this form factor tension on neutrino-
nucleus cross section predictions is studied using nuclear many-body
calculations with the GFMC and SF methods in Sec. IV below. The LQCD results
of Refs. [30, 34] lead to nearly indistinguishable cross-section results that
will be denoted “LQCD Bali et al./Park et al.” or “LQCD” below and used for
comparison with the deuterium bubble-chamber analysis of Ref. [65], denoted
“D2 Meyer et al.” or “D2” below.

Figure 3: The fractional contribution of the axial form factor to the one-body
cross section for MiniBooNE kinematics in selected bins of $\cos\theta_{\mu}$
predicted using the SF formalism shown in Fig. 4. The central value and
uncertainty of the axial form factor correspond to the D2 Meyer et al. $z$
expansion results [65]. Figures are cut off at $T_{\mu}=1$ GeV because the
cross sections become extremely small beyond this point and relative cross-
section uncertainties become large.

##  IV Flux-Averaged Cross Section Results

To evaluate both the nuclear model and nucleon axial form factor dependence of
neutrino-nucleus cross-section predictions and their agreement with data, the
GFMC and spectral function methods are used to predict flux-averaged cross
sections that can be compared with data from the T2K and MiniBooNE
experiments. The MiniBooNE data for this comparison is a double differential
CCQE measurement where the main CC1$\pi$\+ background has been subtracted
using a tuned model [13], and the T2K data is a double differential CC0$\pi$
measurement [114]. Muon neutrino flux-averaged cross sections were calculated
from

$\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}=\int
dE_{\nu}\phi(E_{\nu})\frac{d\sigma(E_{\nu})}{dT_{\mu}d\cos\theta_{\mu}},$ (43)

where $\phi(E_{\nu})$ are the normalized $\nu_{\mu}$ fluxes from MiniBooNE and
T2K. Details on the neutrino fluxes for each experiment can be found in the
references above. $\frac{d\sigma(E_{\nu})}{dT_{\mu}d\cos\theta_{\mu}}$ are the
corresponding inclusive cross sections computed using the GFMC and SF methods
as described in Sec. II.

The fractional contribution of the axial form factor to the one-body piece of
the MiniBooNE flux-averaged cross section is determined by including only pure
axial and axial-vector interference terms in the cross section and shown in
Fig. 3. These pure axial and axial-vector interference terms account for half
or more of the total one-body cross section for most $T_{\mu}$ and
$\cos\theta_{\mu}$, which emphasizes the need for an accurate determination of
the nucleon axial form factor.

Figures 4 and 5 show the GFMC and SF predictions for MiniBooNE and T2K,
respectively, including the breakdown into one-body and two-body
contributions. For these comparisons we use the D2 Meyer et al. $z$ expansion
for $F_{A}$. Two features of the calculations should be noted before
discussing the results of these comparisons. First, the uncertainty bands in
the SF come only from the axial form factor, while the GFMC error bands
include axial form factor uncertainties as well as a combination of GFMC
statistical errors and uncertainties associated with the maximum-entropy
inversion. Secondly, the axial form factor enters into the SF only in the one-
body term, in contrast to the GFMC prediction where it enters into both the
one-body and one and two-body interference term.

Below in Table 1 we quantify the differences between GFMC and SF predictions
for both MiniBooNE and T2K. The percent difference in the differential cross
sections at each model’s peak are shown. The GFMC predictions are up to $20\%$
larger in backwards angle regions for MiniBooNE and $13\%$ larger for T2K in
the same backward region. The agreement between GFMC and SF predictions is
better at more forward angles but a $5$-$10\%$ difference persists.

Figure 4 shows relatively good agreement between SF predictions and MiniBooNE
cross-section results for a range of kinematics, with SF underestimating data
in the lowest bins of $\cos\theta_{\mu}$ by up to $20\%$ at the peak. The
location of the peak is correctly predicted, and the agreement increases at
forward angles. GFMC predictions for MiniBooNE see a similar deficit of events
at low $\cos\theta_{\mu}$. The enhanced disagreement for GFMC at small
$\cos\theta_{\mu}$ and large $T_{\mu}$ is expected as these bins occur at
larger $Q^{2}$ where relativistic effects become important. Comparisons with
T2K for both SF and GFMC in Fig. 5 are in excellent agreement with the data
even for backwards angles. The GFMC prediction here is expected to be more
accurate as the T2K $\nu_{\mu}$ flux has a smaller tail at high neutrino
energy. It is interesting to note the shift in the location of the two-body
peak between GFMC and SF with both MiniBooNE and T2K fluxes, the GFMC two-body
contribution peaking in roughly the same location as the one-body, whereas the
SF two-body peak is shifted to lower $T_{\mu}$. This is again because the GFMC
two-body curve is dominated by the one- and two-body current interference
which roughly peaks in the same location as the pure one-body term. In the SF
calculation, the one- and two-body current interference term yielding two
nucleon emission has been studied in Refs. [84, 83] and found to be
negligible. Therefore, in this work only the pure two-body contribution has
been included which requires large energy transfer and thus peaks at lower
$T_{\mu}$. Accounting for one- and two-body current interference with a single
nucleon emitted in the final state involves some developments in the SF
formalism and it will be the subject of a future work.

MiniBooNE | $0.2<\cos\theta_{\mu}<0.3$ | $0.5<\cos\theta_{\mu}<0.6$ | $0.8<\cos\theta_{\mu}<0.9$
---|---|---|---
GFMC/SF difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $22.8$ | $20.3$ | $5.6$

T2K | $0.0<\cos\theta_{\mu}<0.6$ | $0.80<\cos\theta_{\mu}<0.85$ | $0.94<\cos\theta_{\mu}<0.98$
---|---|---|---
GFMC/SF difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $13.4$ | $7.3$ | $10.0$

Table 1: Difference in value of
$\frac{d\sigma(E_{\nu})}{dT_{\mu}d\cos\theta_{\mu}}$ at the quasielastic peak
computed using GFMC and SF methods for MiniBooNE and T2K flux-averaged double-
differential cross sections.  Figure 4: Breakdown into one- and two-body
current contributions of the $\nu_{\mu}$ flux-averaged differential cross
sections for MiniBooNE: 1b and 2b denotes one- and two-body current
contributions while 12b denotes the total sum of these contributions. The top
panel shows Spectral Function predictions in three bins of $\cos\theta_{\mu}$
with the one-body contributions in orange, two-body contributions in red, and
the total in blue. The lower panel shows GFMC predictions with the same
breakdown between one- and two-body current contributions, although the two-
body results include interference effects only in the GFMC case. The D2 Meyer
et al. $z$ expansion results for $F_{A}$ are used in both cases [65].  Figure
5: Breakdown into one- and two-body current contributions of the $\nu_{\mu}$
flux-averaged differential cross sections for T2K. The color code is as in
Fig. 4. Figure 6: Percent change in the value of the MiniBooNE flux-averaged
cross section for $0.5<\cos\theta_{\mu}<0.6$ vs. percent change in two
parameters describing $\Delta$ resonance production and decay entering
calculations of two-body current (MEC) effects: $C^{A}_{5}(Q^{2})$ is the
dominant $N\rightarrow\Delta$ transition form factor, and $\Lambda_{R}$
renormalizes the self energy of the $\Delta$ as described in Sec. II.2.

MiniBooNE | $0.2<\cos\theta_{\mu}<0.3$ | $0.5<\cos\theta_{\mu}<0.6$ | $0.8<\cos\theta_{\mu}<0.9$
---|---|---|---
SF Difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $16.3$ | $17.1$ | $9.3$
GFMC Difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $18.6$ | $17.1$ | $12.2$

T2K | $0.0<\cos\theta_{\mu}<0.6$ | $0.80<\cos\theta_{\mu}<0.85$ | $0.94<\cos\theta_{\mu}<0.98$
---|---|---|---
SF difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $15.3$ | $8.2$ | $3.3$
GFMC difference in $d\sigma_{\rm{peak}}$ $(\%)$ | $15.8$ | $8.0$ | $4.6$

Table 2: Percent increase in $\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}$ at
the quasielastic peak between predictions using LQCD Bali et al./Park et al.
$z$ expansion versus D2 Meyer et al. $z$ expansion nucleon axial form factor
results.  Figure 7: The $\nu_{\mu}$ flux-averaged differential cross sections
for MiniBooNE. The top panel shows Spectral Function predictions in three bins
of $\cos\theta_{\mu}$ with the D2 Meyer et al. $z$ expansion $F_{A}$ in blue,
as well as the LQCD Bali et al./Park et al. $z$ expansion $F_{A}$ in green.
The dipole parameterization with $M_{A}=1.0$ GeV is shown without
uncertainties as a black line. The lower panel shows GFMC predictions using
the same set of axial form factors, although in the GFMC case systematic
uncertainties including those arising from inversion of the Euclidean response
functions are included in all results and the $M_{A}=1.0$ GeV dipole form
factor results are therefore shown as a black band. Figure 8: The $\nu_{\mu}$
flux-averaged differential cross sections for T2K. Details are as in Fig. 7.
Figure 9: Percent change in peak value of MiniBooNE flux-averaged cross
section for $0.5<\cos\theta_{\mu}<0.6$ vs. percent change in the z expansion
parameters $a_{k}$. Results are shown for predictions using SF (black) and
GFMC (blue) methods, including the slopes extracted from linear fits. Figure
10:  The ratio between relative uncertainties in MiniBooNE cross sections and
nucleon axial form factor $z$ expansion coefficients, $D_{k}(\theta_{\mu})$,
vs muon scattering angle $\theta_{\mu}$ for $z$ expansion coefficients $a_{k}$
with $k\in\\{0,1,2\\}$. Results for $D_{k}(\theta_{\mu})$ are calculated using
the GFMC method with $T_{\mu}$ set equal to the peak value of the flux-
averaged differential cross section for each $\theta_{\mu}$ and correspond to
the slopes shown for particular $\cos\theta_{\mu}$ bins in Fig. 9.

The two-body current contributions to flux-averaged cross sections are among
the most unconstrained due to a lack of knowledge of the $N\rightarrow\Delta$
transition form factors and other nucleon resonance properties and two-nucleon
correlations that enter the theoretical description of these contributions
discussed in Sec. II. Empirical MEC models [115] and phenomenological models
such as the Valencia model [116] are used to describe two-body current effects
in the neutrino event generators used for oscillation analyses at NOvA,
MINERvA, and T2K, and neutrino experiments must resort to further tuning that
significantly changes the shape and normalization of this 2p2h component in
order to achieve satisfactory agreement between theory and data for flux-
averaged cross sections [59, 117]. The degree to which any realistic model can
reproduce the effect of these experimental tunes remains to be seen. To
investigate the model dependence of two-body current effects, we have chosen
two parameters from the SF two-body current calculation and varied them to
investigate their effect on flux-averaged cross sections. These two parameters
are $C^{A}_{5}(Q^{2}=0)$, which is the dominant coupling444Note that although
$C_{6}(Q^{2})>C_{5}(Q^{2})$ for $Q^{2}\lesssim 1$ GeV2 [118], contributions to
neutrino-nucleus cross sections from $C_{6}$ are suppressed by lepton masses
and therefore sub-dominant. A relation between $C_{6}$ and $C_{5}$ analogous
to Eq. (8) is also predicted by leading order chiral perturbation theory. See
Refs. [119, 85] for more details. appearing in the $N\rightarrow\Delta$ vertex
according to phenomenological models [120, 119, 88, 101, 121, 85, 122] and
pioneering LQCD calculations [123, 118], and $\Lambda_{R}$, which is a
parameter that renormalizes the self energy of the $\Delta$. These parameters
have been chosen because they affect the $\Delta$ piece of the two-body
current, which we have seen provides the largest contribution, as well as
because they are highly unconstrained.

Each parameter was varied by $\pm 5,10\%$ and the effect on the flux-averaged
cross section at the peak of the two-body contribution was computed. The
effect can be seen in Fig. 6 where we have plotted the percent change in the
MiniBooNE cross section versus the percent change in each parameter for
$0.5<\cos\theta_{\mu}<0.6$, $T_{\mu}=325$ MeV. This was fit to a line so that
as in Sec. III.1 the extracted slope is an estimate of the derivative of the
cross section with respect to each parameter. The derivative with respect to
$C^{A}_{5}(0)$ is estimated to be 0.31, meaning that achieving a given cross-
section uncertainty requires $C^{A}_{5}(0)$ to be known with $\lesssim 3$
times that uncertainty. A similar though slightly smaller slope of 0.29 is
found for $\Lambda_{R}$. Current extractions of $C_{5}(0)$ rely on single pion
production data from deuterium bubble chamber experiments [10, 11, 12], and
due to limited statistics model assumptions on the relations between
$N\to\Delta$ transition form factors are typically included to reduce the
number of fit parameters. Depending on the model assumptions used, the
resulting uncertainty on $C_{5}(0)$ is estimated to be 10-15$\%$ in the
analysis of Ref. [122], with similar though slightly less conservative
uncertainties estimated in Refs. [85, 121]. Note that all of these analysis
assume a dipole parameterization of $F_{A}$ as well as modified dipole
parameterizations of $C_{5}^{A}$, and therefore it is possible that these
uncertainties are still underestimated. Even less is known about the
uncertainty in determining $\Lambda_{R}$ [89]. A $15\%$ variation in either
$C_{5}^{A}(0)$ or $\Lambda_{R}$ changes the flux-averaged cross section by
roughly $5\%$, and it will therefore be important to obtain more information
on these parameters in order to achieve few-percent precision on cross-section
predictions.

Focusing now on $F_{A}$, Figs. 7 and 8 compare flux-averaged cross sections
with different axial form factor determinations: a dipole form factor with
$M_{A}=1.0$ GeV, the D2 Meyer et al. $z$ expansion, and the LQCD Bali et
al./Park et al. $z$ expansion. One can see that the LQCD $z$ expansion
increases the normalization of the cross section across the whole phase space,
with significantly more enhancement in the bins of low $\cos\theta_{\mu}$
corresponding to backward angles and higher $Q^{2}$. This is quantified in
Table 2, which shows the percentage difference in the peak values of
$\frac{d\sigma}{dT_{\mu}d\cos\theta_{\mu}}$ for the LQCD and D2 $z$ expansion
results. The LQCD prediction increases the peak cross section between
$10$-$20\%$, with the discrepancy growing at backwards angles.

To investigate the sensitivity of the flux-averaged differential cross section
to variations in the axial form factor, derivatives of the MiniBooNE cross
section with respect to the model-independent $z$ expansion parameters $a_{k}$
are computed as described in Sec. III.1. Figure 9 shows the percent
differences in flux-averaged cross sections evaluated at the quasielastic peak
that have been computed using both GFMC and SF methods after independently
varying each $a_{k}$ by $\pm 5,10\%$. The slopes of the resulting linear fits
provide model-independent determinations of the sensitivity of the peak cross
section to variations in $F_{A}$. It is clear that the impact of varying each
$a_{k}$ decreases as $k$ increases, as expected since the contribution of each
$a_{k}$ is suppressed by the $k$-th power of $z(Q^{2})<1$. In particular, a
$10\%$ change in $a_{0}$ results in a $10\%$ change to the peak cross section,
while a $10\%$ change in $a_{1}$ results in a $1\%$ change in the peak cross
section, and $10\%$ variation of $a_{k}$ with $k\geq 2$ leads to sub-percent
changes in the peak cross section. It is noteworthy that the uncertainty
relations determined using GFMC and SF methods are reasonably consistent, with
slightly larger slopes arising in GFMC calculations where the axial form
factor appears in one- and two-body interference terms in addition to one-body
contributions. We have further studied the $\cos\theta_{\mu}$ dependence of
these uncertainty relations using the GFMC method as shown in Fig. 10. It can
be seen that the uncertainty relation factor $D_{k}(\theta_{\mu})$ defined in
Eq. (41), which is equivalent to the slopes shown for particular
$\cos\theta_{\mu}$ bins in Fig. 9, decreases as $\cos\theta_{\mu}$ increases,
except for $\cos\theta_{\mu}\geq 0.75$ for $a_{1}$ and $a_{2}$. These
inflection points are merely due to the absolute values appearing in the
definition of $D_{k}(\theta_{\mu})$. Achieving $1\%$ precision in $a_{0}$ and
$10\%$ precision in $a_{1}$ is seen to be sufficient to ensure that axial form
factor uncertainties lead to only $1\%$ uncertainties in MiniBooNE flux-
averaged differential cross sections across the entire phase space. The D2
Meyer et al. results lead to $6\%$ precision on $a_{1}$ but only $10\%$
precision on $a_{0}$,555Note that $a_{0}$ is highly correlated with the other
$a_{k}$ in this determination because $a_{0}$ is fixed by treating Eq. (36) as
a constraint. and therefore higher-precision determinations of the axial form
factor are necessary for achieving few-percent cross-section uncertainties.
Current LQCD results quote higher precision on $a_{0}$ [34, 35], though the
LQCD studies use different values of $t_{0}$ that prevent an immediate
quantitative comparison with these uncertainty targets. It would simplify
comparisons and uncertainty quantification for future LQCD calculations to use
the Meyer et al. value of $t_{0}=-0.28$ GeV2, which is chosen to optimize the
rate of convergence of the $z$ expansion over the interval $0\leq Q^{2}\leq 1$
GeV2 [65].

##  V Conclusions

Precise theoretical calculations of neutrino scattering cross sections on
target nuclei supplemented by reliable estimates of the associated
uncertainties are essential for the success of the accelerator neutrino
experimental program. A promising strategy for achieving this goal is to
combine state-of-the-art nuclear many-body methods with LQCD calculations of
few-nucleon observables, including single-nucleon form factors. In addition,
it is critical to quantify the precision level required on such input
parameters to compute neutrino-nucleus cross section with a given accuracy.

In this work, we analyze the impact of varying the nucleon axial form factors
on flux-averaged charged-current neutrino-nucleus cross-sections as measured
by the MiniBooNE and T2K experiments. Specifically, we consider a dipole form
factor with $M_{A}=1.0$ GeV, the D2 Meyer et al. $z$ expansion results fit to
experimental deuterium bubble-chamber neutrinos scattering data [65], and the
LQCD Bali et al./Park et al. $z$ expansion results [30, 34], which are
consistent with other recent LQCD calculations [35]. The nuclear many-body
problem is solved using the GFMC and SF approaches, which rely on similar
descriptions of the initial target state. For instance, the one- and two-
nucleon SF of 12C are obtained from state-of-the-art VMC calculations that use
as input the same nuclear Hamiltonian as the GFMC.

For the first time, we validate the sum of one- and two-body current
contributions obtained within the SF approach against flux-averaged
electroweak cross sections. Using the D2 Meyer et al. form factors, we observe
an overall good agreement with the MiniBooNE and T2K data — there are tensions
with the MiniBooNE data for small values of $\cos\theta_{\mu}$. At small
$T_{\mu}$, these discrepancies might be due to the model dependent procedure
adopted to subtract the pion-absorption contribution from the experimental
data [13]. Using instead the LQCD axial form factor in the SF calculations
show an improved agreement with MiniBooNE data for these kinematics over the
entire $T_{\mu}$ range. By contrast, the full SF predictions (with one- and
two-body currents) using either choice of axial form factors appear to provide
a good description of the T2K data for all the kinematics that we consider.
Discrepancies between GFMC results and MiniBooNE data are seen for small
values of $\cos\theta_{\mu}$, independent of the form factor of choice. They
could be due to the nonrelativistic nature of the GFMC calculation in
conjunction with the significant large-$E_{\nu}$ tail of the MiniBooNE flux.
The GFMC results similarly show good agreement with T2K data using both form
factor parameterizations. It is noteworthy that adopting the LQCD $z$
expansion results leads to a $10$-$20\%$ enhancement of the GFMC and SF cross
section with respect to the D2 $z$ expansion across the whole phase space,
with the largest effect seen in bins of low $\cos\theta_{\mu}$. We find that
the SF and GFMC methods predict similar uncertainty-quantification results on
the form factor dependence. The sensitivity analysis that we carry out in this
work enables to determine an uncertainty targets for future LQCD axial form
factor studies.

We observe that the GFMC and SF cross sections including one- and two-body
contributions are consistent for the different kinematics and data sets
considered. The differences in the strength and peak positions arise from
different factors such as the nonrelativistic nature of the GFMC calculations
and the static treatment of the $\Delta$ propagator in the nonrelativistic
reduction of the two-body currents. The use of a factorization scheme that
neglects final state interactions between the emitted particles and the
spectator nucleons and the lack of one- and two-body current interference in
the SF approach also contribute to the observed discrepancies. The inclusion
of relativistic corrections in the kinematics of the GFMC electroweak
responses is currently ongoing and it will be the subject of a future work;
modeling one- and two-body interference effects yielding one nucleon emission
in the final state requires some conceptual developments in the SF approach
that will be investigated in the future.

Reaction channels yielding the production of one or more pions in the final
state have not been considered here. While including explicit pion degrees of
freedom in the GFMC involves nontrivial difficulties [124], both resonant and
non-resonant pion production mechanisms have been included in the SF method
and validated against electron scattering data [24]. We plan to carry out an
analysis of the cross-section dependence on the input parameters needed to
describe non-resonant $N\to N\pi$ transition amplitudes in a future work.
Pioneering LQCD calculations of $N\to N\pi$ transition amplitudes have begun
[39] and will eventually supply robust QCD predictions to constrain these
input parameters. It is also noteworthy that there has been significant
efforts to calculate nucleon properties using quark models involving
approximate solutions to QCD Schwinger-Dyson and Fadeev equations [125, 126]
that can reproduce LQCD results for nucleon axial and pseudoscalar form
factors [56, 57] and could provide valuable predictions of resonant and non-
resonant pion production amplitudes that are not yet known precisely from
LQCD.

In order to provide predictions for 16O and 40Ar needed for next-generation
accelerator neutrino experiments, we will employ QMC methods that are capable
to treat larger nuclear systems, e.g., the auxiliary field Diffusion Monte
Carlo method [127] and assess the uncertainty of our predictions. We plan to
validate future predictions for 40Ar against accurate data from the Short-
Baseline Neutrino Program as they become available and continue studying the
nuclear many-body and nucleon form factor precision required to achieve the
design cross-section uncertainty goals of DUNE.

###### Acknowledgements.

We thank Gunnar Bali, Rajan Gupta, Andreas Kronfeld, and Aaron Meyer for
helpful discussions. This manuscript has been authored by Fermi Research
Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of
Energy, Office of Science, Office of High Energy Physics and Fermilab LDRD
awards (N.S). The present research is also supported by the U.S. Department of
Energy, Office of Science, Office of Nuclear Physics, under contracts DE-
AC02-06CH11357 (A.L.), by the NUCLEI SciDAC program (A.L.), the DOE Early
Career Research Program (A.L.), and by the Visiting Scholars Award Program of
the Universities Research Association #21-S-12. Quantum Monte Carlo
calculations were performed on the parallel computers of the Laboratory
Computing Resource Center, Argonne National Laboratory, the computers of the
Argonne Leadership Computing Facility via the INCITE grant “Ab-initio nuclear
structure and nuclear reactions”.

## References

  * Acciarri  _et al._ [2015] R. Acciarri _et al._ (DUNE), (2015), arXiv:1512.06148 [physics.ins-det] .
  * Abe  _et al._ [2021] K. Abe _et al._ (T2K), Phys. Rev. D 103, 112008 (2021), arXiv:2101.03779 [hep-ex] .
  * Acero  _et al._ [2018] M. A. Acero _et al._ (NOvA), Phys. Rev. D 98, 032012 (2018), arXiv:1806.00096 [hep-ex] .
  * Mann  _et al._ [1973] W. A. Mann _et al._ , Phys. Rev. Lett. 31, 844 (1973).
  * Barish  _et al._ [1977] S. J. Barish _et al._ , Phys. Rev. D 16, 3103 (1977).
  * Baker  _et al._ [1981] N. J. Baker, A. M. Cnops, P. L. Connolly, S. A. Kahn, H. G. Kirk, M. J. Murtagh, R. B. Palmer, N. P. Samios, and M. Tanaka, Phys. Rev. D 23, 2499 (1981).
  * Miller  _et al._ [1982] K. L. Miller _et al._ , Phys. Rev. D 26, 537 (1982).
  * Kitagaki  _et al._ [1983] T. Kitagaki _et al._ , Phys. Rev. D 28, 436 (1983).
  * Barish  _et al._ [1979] S. J. Barish _et al._ , Phys. Rev. D 19, 2521 (1979).
  * Radecky  _et al._ [1982] G. M. Radecky _et al._ , Phys. Rev. D 25, 1161 (1982), [Erratum: Phys.Rev.D 26, 3297 (1982)].
  * Kitagaki  _et al._ [1986] T. Kitagaki _et al._ , Phys. Rev. D 34, 2554 (1986).
  * Kitagaki  _et al._ [1990] T. Kitagaki _et al._ , Phys. Rev. D 42, 1331 (1990).
  * Aguilar-Arevalo  _et al._ [2010] A. A. Aguilar-Arevalo _et al._ (MiniBooNE), Phys. Rev. D 81, 092005 (2010), arXiv:1002.2680 [hep-ex] .
  * Benhar and Meloni [2009] O. Benhar and D. Meloni, Phys. Rev. D 80, 073003 (2009), arXiv:0903.2329 [hep-ph] .
  * Martini  _et al._ [2010] M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 81, 045502 (2010), arXiv:1002.4538 [hep-ph] .
  * Martini  _et al._ [2011] M. Martini, M. Ericson, and G. Chanfray, Phys. Rev. C 84, 055502 (2011), arXiv:1110.0221 [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] .
  * Nieves  _et al._ [2012] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B 707, 72 (2012), arXiv:1106.5374 [hep-ph] .
  * Megias  _et al._ [2015] G. D. Megias _et al._ , Phys. Rev. D 91, 073004 (2015), arXiv:1412.1822 [nucl-th] .
  * Gonzaléz-Jiménez  _et al._ [2014] R. Gonzaléz-Jiménez, G. D. Megias, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. C 90, 035501 (2014), arXiv:1407.8346 [nucl-th] .
  * Pandey  _et al._ [2015] V. Pandey, N. Jachowicz, T. Van Cuyck, J. Ryckebusch, and M. Martini, Phys. Rev. C 92, 024606 (2015), arXiv:1412.4624 [nucl-th] .
  * Lovato  _et al._ [2020] A. Lovato, J. Carlson, S. Gandolfi, N. Rocco, and R. Schiavilla, Phys. Rev. X 10, 031068 (2020), arXiv:2003.07710 [nucl-th] .
  * Rocco  _et al._ [2019a] N. Rocco, C. Barbieri, O. Benhar, A. De Pace, and A. Lovato, Phys. Rev. C 99, 025502 (2019a), arXiv:1810.07647 [nucl-th] .
  * Rocco  _et al._ [2019b] N. Rocco, S. X. Nakamura, T. S. H. Lee, and A. Lovato, Phys. Rev. C 100, 045503 (2019b), arXiv:1907.01093 [nucl-th] .
  * Andreoli  _et al._ [2022] L. Andreoli, J. Carlson, A. Lovato, S. Pastore, N. Rocco, and R. B. Wiringa, Phys. Rev. C 105, 014002 (2022), arXiv:2108.10824 [nucl-th] .
  * Korover  _et al._ [2022] I. Korover _et al._ (CLAS), (2022), arXiv:2209.01492 [nucl-ex] .
  * Shintani  _et al._ [2019] E. Shintani, K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, and T. Yamazaki, Phys. Rev. D 99, 014510 (2019), [Erratum: Phys.Rev.D 102, 019902 (2020)], arXiv:1811.07292 [hep-lat] .
  * Jang  _et al._ [2020a] Y.-C. Jang, R. Gupta, B. Yoon, and T. Bhattacharya, Phys. Rev. Lett. 124, 072002 (2020a), arXiv:1905.06470 [hep-lat] .
  * Jang  _et al._ [2020b] Y.-C. Jang, R. Gupta, H.-W. Lin, B. Yoon, and T. Bhattacharya, Phys. Rev. D 101, 014507 (2020b), arXiv:1906.07217 [hep-lat] .
  * Bali  _et al._ [2020] G. S. Bali, L. Barca, S. Collins, M. Gruber, M. Löffler, A. Schäfer, W. Söldner, P. Wein, S. Weishäupl, and T. Wurm (RQCD), JHEP 05, 126, arXiv:1911.13150 [hep-lat] .
  * Alexandrou  _et al._ [2021] C. Alexandrou _et al._ , Phys. Rev. D 103, 034509 (2021), arXiv:2011.13342 [hep-lat] .
  * Ishikawa  _et al._ [2021] K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, E. Shintani, and T. Yamazaki (PACS), Phys. Rev. D 104, 074514 (2021), arXiv:2107.07085 [hep-lat] .
  * Meyer  _et al._ [2022a] A. S. Meyer _et al._ , PoS LATTICE2021, 081 (2022a), arXiv:2111.06333 [hep-lat] .
  * Park  _et al._ [2022] S. Park, R. Gupta, B. Yoon, S. Mondal, T. Bhattacharya, Y.-C. Jang, B. Joó, and F. Winter (Nucleon Matrix Elements (NME)), Phys. Rev. D 105, 054505 (2022), arXiv:2103.05599 [hep-lat] .
  * Djukanovic  _et al._ [2022] D. Djukanovic, G. von Hippel, J. Koponen, H. B. Meyer, K. Ottnad, T. Schulz, and H. Wittig, (2022), arXiv:2207.03440 [hep-lat] .
  * Alexandrou  _et al._ [2019] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen, G. Koutsou, and A. Vaquero Aviles-Casco, Phys. Rev. D 100, 014509 (2019), arXiv:1812.10311 [hep-lat] .
  * Bar [2019a] O. Bar, Phys. Rev. D 99, 054506 (2019a), arXiv:1812.09191 [hep-lat] .
  * Bar [2019b] O. Bar, Phys. Rev. D 100, 054507 (2019b), arXiv:1906.03652 [hep-lat] .
  * Barca  _et al._ [2022] L. Barca, G. S. Bali, and S. Collins, PoS LATTICE2021, 359 (2022), arXiv:2110.11908 [hep-lat] .
  * Peardon  _et al._ [2009] M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Joo, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), Phys. Rev. D 80, 054506 (2009), arXiv:0905.2160 [hep-lat] .
  * Morningstar  _et al._ [2011] C. Morningstar, J. Bulava, J. Foley, K. J. Juge, D. Lenkner, M. Peardon, and C. H. Wong, Phys. Rev. D 83, 114505 (2011), arXiv:1104.3870 [hep-lat] .
  * Dudek  _et al._ [2013a] J. J. Dudek, R. G. Edwards, and C. E. Thomas (Hadron Spectrum), Phys. Rev. D 87, 034505 (2013a), [Erratum: Phys.Rev.D 90, 099902 (2014)], arXiv:1212.0830 [hep-ph] .
  * Dudek  _et al._ [2013b] J. J. Dudek, R. G. Edwards, P. Guo, and C. E. Thomas (Hadron Spectrum), Phys. Rev. D 88, 094505 (2013b), arXiv:1309.2608 [hep-lat] .
  * Wilson  _et al._ [2015] D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D 92, 094502 (2015), arXiv:1507.02599 [hep-ph] .
  * Lang and Verduci [2013] C. B. Lang and V. Verduci, Phys. Rev. D 87, 054502 (2013), arXiv:1212.5055 [hep-lat] .
  * Kiratidis  _et al._ [2015] A. L. Kiratidis, W. Kamleh, D. B. Leinweber, and B. J. Owen, Phys. Rev. D 91, 094509 (2015), arXiv:1501.07667 [hep-lat] .
  * Andersen  _et al._ [2018] C. W. Andersen, J. Bulava, B. Hörz, and C. Morningstar, Phys. Rev. D 97, 014506 (2018), arXiv:1710.01557 [hep-lat] .
  * Silvi  _et al._ [2021] G. Silvi _et al._ , Phys. Rev. D 103, 094508 (2021), arXiv:2101.00689 [hep-lat] .
  * Bulava  _et al._ [2022] J. Bulava, A. D. Hanlon, B. Hörz, C. Morningstar, A. Nicholson, F. Romero-López, S. Skinner, P. Vranas, and A. Walker-Loud, (2022), arXiv:2208.03867 [hep-lat] .
  * Francis  _et al._ [2019] A. Francis, J. R. Green, P. M. Junnarkar, C. Miao, T. D. Rae, and H. Wittig, Phys. Rev. D 99, 074505 (2019), arXiv:1805.03966 [hep-lat] .
  * Hörz  _et al._ [2021] B. Hörz _et al._ , Phys. Rev. C 103, 014003 (2021), arXiv:2009.11825 [hep-lat] .
  * Amarasinghe  _et al._ [2021] S. Amarasinghe, R. Baghdadi, Z. Davoudi, W. Detmold, M. Illa, A. Parreno, A. V. Pochinsky, P. E. Shanahan, and M. L. Wagman, (2021), arXiv:2108.10835 [hep-lat] .
  * Kronfeld  _et al._ [2019] A. S. Kronfeld, D. G. Richards, W. Detmold, R. Gupta, H.-W. Lin, K.-F. Liu, A. S. Meyer, R. Sufian, and S. Syritsyn (USQCD), Eur. Phys. J. A 55, 196 (2019), arXiv:1904.09931 [hep-lat] .
  * Meyer  _et al._ [2022b] A. S. Meyer, A. Walker-Loud, and C. Wilkinson 10.1146/annurev-nucl-010622-120608 (2022b), arXiv:2201.01839 [hep-lat] .
  * Ruso  _et al._ [2022] L. A. Ruso _et al._ , (2022), arXiv:2203.09030 [hep-ph] .
  * Chen  _et al._ [2022] C. Chen, C. S. Fischer, C. D. Roberts, and J. Segovia, Phys. Rev. D 105, 094022 (2022), arXiv:2103.02054 [hep-ph] .
  * Chen and Roberts [2022] C. Chen and C. D. Roberts, (2022), arXiv:2206.12518 [hep-ph] .
  * Andreopoulos  _et al._ [2010] C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Acero  _et al._ [2020] M. A. Acero _et al._ (NOvA), Eur. Phys. J. C 80, 1119 (2020), arXiv:2006.08727 [hep-ex] .
  * Tena-Vidal  _et al._ [2021] J. Tena-Vidal _et al._ (GENIE), Phys. Rev. D 104, 072009 (2021), arXiv:2104.09179 [hep-ph] .
  * Tena-Vidal  _et al._ [2022] J. Tena-Vidal _et al._ (GENIE), (2022), arXiv:2206.11050 [hep-ph] .
  * Hill [2006] R. J. Hill, eConf C060409, 027 (2006), arXiv:hep-ph/0606023 .
  * Hill and Paz [2010] R. J. Hill and G. Paz, Phys. Rev. D 82, 113005 (2010), arXiv:1008.4619 [hep-ph] .
  * 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] .
  * Nakamura  _et al._ [2002] S. Nakamura, T. Sato, S. Ando, T. S. Park, F. Myhrer, V. P. Gudkov, and K. Kubodera, Nucl. Phys. A 707, 561 (2002), arXiv:nucl-th/0201062 .
  * Nakamura  _et al._ [2010] K. Nakamura _et al._ (Particle Data Group), J. Phys. G 37, 075021 (2010).
  * Shen  _et al._ [2012] G. Shen, L. E. Marcucci, J. Carlson, S. Gandolfi, and R. Schiavilla, Phys. Rev. C 86, 035503 (2012), arXiv:1205.4337 [nucl-th] .
  * Wiringa  _et al._ [1995] R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. Rev. C 51, 38 (1995), arXiv:nucl-th/9408016 .
  * Pieper [2008] S. C. Pieper, AIP Conf. Proc. 1011, 143 (2008).
  * Bryan [1990] R. Bryan, Eur. Biophys. J. 18, 165 (1990).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. E. Gubernatis, Phys. Rept. 269, 133 (1996).
  * Lovato  _et al._ [2016] A. Lovato, S. Gandolfi, J. Carlson, S. C. Pieper, and R. Schiavilla, Phys. Rev. Lett. 117, 082501 (2016), arXiv:1605.00248 [nucl-th] .
  * Lovato  _et al._ [2018] A. Lovato, S. Gandolfi, J. Carlson, E. Lusk, S. C. Pieper, and R. Schiavilla, Phys. Rev. C 97, 022502 (2018), arXiv:1711.02047 [nucl-th] .
  * Raghavan  _et al._ [2021] K. Raghavan, P. Balaprakash, A. Lovato, N. Rocco, and S. M. Wild, Phys. Rev. C 103, 035502 (2021), arXiv:2010.12703 [nucl-th] .
  * Rocco  _et al._ [2018] N. Rocco, W. Leidemann, A. Lovato, and G. Orlandini, Phys. Rev. C 97, 055501 (2018), arXiv:1801.07111 [nucl-th] .
  * Barrow  _et al._ [2021] J. L. Barrow, S. Gardiner, S. Pastore, M. Betancourt, and J. Carlson, Phys. Rev. D 103, 052001 (2021), arXiv:2010.04154 [nucl-th] .
  * Bernard  _et al._ [1994] V. Bernard, N. Kaiser, and U.-G. Meissner, Phys. Rev. D 50, 6899 (1994).
  * Bernard  _et al._ [2002] V. Bernard, L. Elouadrhiri, and U.-G. Meissner, J. Phys. G 28, R1 (2002), arXiv:hep-ph/0107088 .
  * Andreev  _et al._ [2007] V. A. Andreev, T. I. Banks, T. A. Case, D. B. Chitwood, S. M. Clayton, K. M. Crowe, J. Deutsch, J. Egger, S. J. Freedman, V. A. Ganzha, T. Gorringe, F. E. Gray, D. W. Hertzog, M. Hildebrandt, P. Kammel, B. Kiburg, S. Knaack, P. A. Kravtsov, A. G. Krivshich, B. Lauss, K. L. Lynch, E. M. Maev, O. E. Maev, F. Mulhauser, C. S. Özben, C. Petitjean, G. E. Petrov, R. Prieels, G. N. Schapkin, G. G. Semenchuk, M. A. Soroka, V. Tishchenko, A. A. Vasilyev, A. A. Vorobyov, M. E. Vznuzdaev, and P. Winter (MuCap Collaboration), Phys. Rev. Lett. 99, 032002 (2007).
  * Carlson and Schiavilla [1998] J. Carlson and R. Schiavilla, Rev. Mod. Phys. 70, 743 (1998).
  * Kelly [2004] J. J. Kelly, Phys. Rev. C 70, 068202 (2004).
  * Benhar  _et al._ [2015] O. Benhar, A. Lovato, and N. Rocco, Phys. Rev. C 92, 024602 (2015), arXiv:1502.00887 [nucl-th] .
  * Rocco  _et al._ [2016] N. Rocco, A. Lovato, and O. Benhar, Phys. Rev. Lett. 116, 192501 (2016), arXiv:1512.07426 [nucl-th] .
  * Hernandez  _et al._ [2007] E. Hernandez, J. Nieves, and M. Valverde, Phys. Rev. D 76, 033005 (2007), arXiv:hep-ph/0701149 .
  * Ruiz Simo  _et al._ [2017] I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero, and T. W. Donnelly, J. Phys. G 44, 065105 (2017), arXiv:1604.08423 [nucl-th] .
  * Lalakulich  _et al._ [2006] O. Lalakulich, E. A. Paschos, and G. Piranishvili, Phys. Rev. D 74, 014009 (2006), arXiv:hep-ph/0602210 .
  * Paschos  _et al._ [2004] E. A. Paschos, J.-Y. Yu, and M. Sakuda, Phys. Rev. D 69, 014013 (2004), arXiv:hep-ph/0308130 .
  * Dekker  _et al._ [1994] M. J. Dekker, P. J. Brussaard, and J. A. Tjon, Phys. Rev. C 49, 2650 (1994).
  * De Pace  _et al._ [2003] A. De Pace, M. Nardi, W. M. Alberico, T. W. Donnelly, and A. Molinari, Nucl. Phys. A 726, 303 (2003), arXiv:nucl-th/0304084 .
  * Lee [1983] T. s. h. Lee, Phys. Rev. Lett. 50, 1571 (1983).
  * Lee [1984] T. S. H. Lee, Phys. Rev. C 29, 195 (1984).
  * Lee and Matsuyama [1985] T. S. H. Lee and A. Matsuyama, Phys. Rev. C 32, 516 (1985).
  * Lee and Matsuyama [1987] T. S. H. Lee and A. Matsuyama, Phys. Rev. C 36, 1459 (1987).
  * [95] VMC and GFMC radial overlaps, https://www.phy.anl.gov/theory/research/overlap_gfmc/index.html.
  * [96] VMC two-body momentum distributions, https://www.phy.anl.gov/theory/research/momenta2/.
  * Wiringa  _et al._ [2014] R. B. Wiringa, R. Schiavilla, S. C. Pieper, and J. Carlson, Phys. Rev. C 89, 024305 (2014), arXiv:1309.3794 [nucl-th] .
  * Weiss  _et al._ [2015] R. Weiss, B. Bazak, and N. Barnea, Phys. Rev. C 92, 054311 (2015), arXiv:1503.07047 [nucl-th] .
  * Weiss  _et al._ [2019] R. Weiss, I. Korover, E. Piasetzky, O. Hen, and N. Barnea, Phys. Lett. B 791, 242 (2019), arXiv:1806.10217 [nucl-th] .
  * Lonardoni  _et al._ [2018] D. Lonardoni, S. Gandolfi, X. B. Wang, and J. Carlson, Phys. Rev. C 98, 014322 (2018), arXiv:1804.08027 [nucl-th] .
  * Leitner  _et al._ [2009] T. Leitner, O. Buss, L. Alvarez-Ruso, and U. Mosel, Phys. Rev. C 79, 034601 (2009), arXiv:0812.0587 [nucl-th] .
  * Dai  _et al._ [2018] H. Dai _et al._ (Jefferson Lab Hall A), Phys. Rev. C 98, 014617 (2018), arXiv:1803.01910 [nucl-ex] .
  * Barbieri  _et al._ [2019] C. Barbieri, N. Rocco, and V. Somà, Phys. Rev. C 100, 062501 (2019), arXiv:1907.01122 [nucl-th] .
  * Ankowski and Friedland [2020] A. M. Ankowski and A. Friedland, Phys. Rev. D 102, 053001 (2020), arXiv:2006.11944 [hep-ph] .
  * Khachatryan  _et al._ [2021] M. Khachatryan _et al._ (CLAS, e4v), Nature 599, 565 (2021).
  * Patrignani  _et al._ [2016] C. Patrignani _et al._ (Particle Data Group), Chin. Phys. C 40, 100001 (2016).
  * Bodek  _et al._ [2008] A. Bodek, S. Avvakumov, R. Bradford, and H. S. Budd, Eur. Phys. J. C 53, 349 (2008), arXiv:0708.1946 [hep-ex] .
  * Lepage and Brodsky [1980] G. P. Lepage and S. J. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * Lee  _et al._ [2015] G. Lee, J. R. Arrington, and R. J. Hill, Phys. Rev. D 92, 013013 (2015), arXiv:1505.01489 [hep-ph] .
  * Singh [1972] S. K. Singh, Nucl. Phys. B 36, 419 (1972).
  * Alvarez-Ruso  _et al._ [2022] L. Alvarez-Ruso _et al._ , (2022), arXiv:2203.11298 [hep-ex] .
  * Liu  _et al._ [1994] K. F. Liu, S. J. Dong, T. Draper, J. M. Wu, and W. Wilcox, Phys. Rev. D 49, 4755 (1994), arXiv:hep-lat/9305025 .
  * Alexandrou  _et al._ [2007a] C. Alexandrou, G. Koutsou, T. Leontiou, J. W. Negele, and A. Tsapalis, Phys. Rev. D 76, 094511 (2007a), [Erratum: Phys.Rev.D 80, 099901 (2009)], arXiv:0706.3011 [hep-lat] .
  * Abe  _et al._ [2016] K. Abe _et al._ (T2K), Phys. Rev. D 93, 112012 (2016), arXiv:1602.03652 [hep-ex] .
  * Katori [2015] T. Katori, AIP Conf. Proc. 1663, 030001 (2015), arXiv:1304.6014 [nucl-th] .
  * Gran  _et al._ [2013] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D 88, 113007 (2013), arXiv:1307.8105 [hep-ph] .
  * Ascencio  _et al._ [2022] M. V. Ascencio _et al._ (MINERvA), Phys. Rev. D 106, 032001 (2022), arXiv:2110.13372 [hep-ex] .
  * Alexandrou  _et al._ [2011] C. Alexandrou, G. Koutsou, J. W. Negele, Y. Proestos, and A. Tsapalis, Phys. Rev. D 83, 014501 (2011), arXiv:1011.3233 [hep-lat] .
  * Llewellyn Smith [1972] C. H. Llewellyn Smith, Phys. Rept. 3, 261 (1972).
  * Adler [1968] S. L. Adler, Annals Phys. 50, 189 (1968).
  * Graczyk  _et al._ [2009] K. M. Graczyk, D. Kielczewska, P. Przewlocki, and J. T. Sobczyk, Phys. Rev. D 80, 093001 (2009), arXiv:0908.2175 [hep-ph] .
  * Hernandez  _et al._ [2010] E. Hernandez, J. Nieves, M. Valverde, and M. J. Vicente Vacas, Phys. Rev. D 81, 085046 (2010), arXiv:1001.4416 [hep-ph] .
  * Alexandrou  _et al._ [2007b] C. Alexandrou, T. Leontiou, J. W. Negele, and A. Tsapalis, Phys. Rev. Lett. 98, 052003 (2007b), arXiv:hep-lat/0607030 .
  * Madeira  _et al._ [2018] L. Madeira, A. Lovato, F. Pederiva, and K. E. Schmidt, Phys. Rev. C 98, 034005 (2018), arXiv:1803.10725 [nucl-th] .
  * Barabanov  _et al._ [2021] M. Y. Barabanov _et al._ , Prog. Part. Nucl. Phys. 116, 103835 (2021), arXiv:2008.07630 [hep-ph] .
  * Brodsky  _et al._ [2020] S. J. Brodsky _et al._ , Int. J. Mod. Phys. E 29, 2030006 (2020), arXiv:2006.06802 [hep-ph] .
  * Schmidt and Fantoni [1999] K. E. Schmidt and S. Fantoni, Phys. Lett. B446, 99 (1999).


# A Note on Fermat’s Last Theorem for $n=4$

Matan Eliashar and Nati Linial

Fermat’s last theorem (FLT) famously states that the equation
$x^{n}+y^{n}=z^{n}$ has no solution in positive integers $x,y,z$ for any
integer exponent $n>2$. But does this theorem have a quantitative version?
Specifically, for a given $n\geq 3$, how large can $d_{n}(\cdot)$ be if it is
true that $|z^{n}-x^{n}-y^{n}|\geq d_{n}(z)$ for every three positive integers
$z>x,y>0$?

Exponent $n=4$ is the easiest case of FLT and was already known to Fermat. In
fact, his proof shows that even the equation $x^{4}+y^{4}=z^{2}$ has no
solution in positive integers $x,y,z$, where $z>x^{2},y^{2}$ [1]. We consider
here approximate solutions of this equation. A computer search has yielded the
following triplets:

n | 0 | 1 | 2 | 3
---|---|---|---|---
x | 22 | 1058 | 50806 | 2439746
y | 23 | 1103 | 52967 | 2543519
z | 717 | 1653213 | 3812308653 | 8791182100413

with

$x_{n}^{4}+y_{n}^{4}-8=z_{n}^{2}$ (1)

These turn out to be the initial terms of an infinite sequence of triplets
$\left(x_{n},y_{n},z_{n}\right)_{n=1}^{\infty}$, that satisfy Equation (1).
Here $x_{n},y_{n}$ are defined by the above initial conditions and the
recurrences

$x_{n}=48\cdot x_{n-1}+x_{n-2}~{}~{};~{}~{}y_{n}=48\cdot y_{n-1}+y_{n-2}$ (2)

whereas

$z_{n}=2306\cdot z_{n-1}-z_{n-2}+(-1)^{n}\cdot 192$ (3)

These recurrences and the initial conditions clearly imply that
$x_{n},y_{n},z_{n}$ are always positive integers. Using the three-term
recurrence (2) and the initial conditions, we get the closed-form expressions
for $x_{n},y_{n}$:

$x_{n}=a\lambda_{1}^{n}+b\lambda_{2}^{n}~{}~{};~{}~{}y_{n}=c\lambda_{1}^{n}+d\lambda_{2}^{n}$
(4)

where

$\lambda_{1}=\left(24+\sqrt{577}\right)~{},~{}\lambda_{2}=\left(24-\sqrt{577}\right)$

are the two roots of the quadratic $\lambda^{2}=48\lambda+1$ corresponding to
the recurrence (2), and

$a=\left(11+\frac{265}{\sqrt{577}}\right)\;,\;b=\left(11-\frac{265}{\sqrt{577}}\right)\;,\;c=\left(\frac{23}{2}+\frac{551}{2\sqrt{577}}\right)\;,\;d=\left(\frac{23}{2}-\frac{551}{2\sqrt{577}}\right)$

Note also that $\lambda_{1}\cdot\lambda_{2}=-1$, the constant term in that
quadratic.
Consequently:

$x_{n}^{4}=a^{4}\lambda_{1}^{4n}+\left(-1\right)^{n}4a^{3}b\lambda_{1}^{2n}+6a^{2}b^{2}+\left(-1\right)^{n}4ab^{3}\lambda_{1}^{-2n}+b^{4}\lambda_{1}^{-4n}$

and hence:

$x_{n}^{4}+y_{n}^{4}-8=$
$\left(a^{4}+c^{4}\right)\lambda_{1}^{4n}+\left(-1\right)^{n}\left(4a^{3}b+4c^{3}d\right)\lambda_{1}^{2n}+\left(6a^{2}b^{2}+6c^{2}d^{2}-8\right)+\left(-1\right)^{n}\left(4ab^{3}+4cd^{3}\right)\lambda_{1}^{-2n}+\left(b^{4}+d^{4}\right)\lambda_{1}^{-4n}$

The recurrence (3) that $z_{n}$ satisfies is a bit less standard, but its
closed-form solution is not hard to derive. Namely,

$z_{n}=e\mu_{1}^{n}+f\mu_{2}^{n}+(-1)^{n}\cdot g$ (5)

where

$\mu_{1}=\left(1153+48\cdot\sqrt{577}\right)~{},~{}\mu_{2}=\left(1153-48\cdot\sqrt{577}\right)$

are the two roots of the quadratic $\mu^{2}=2306\mu-1$ corresponding to three-
term part of recurrence (3), and

$e=\left(\frac{413661}{1154}+\frac{17221}{1154}\sqrt{577}\right)\;,\;f=\left(\frac{413661}{1154}-\frac{17221}{1154}\sqrt{577}\right)\;,\;g=\left(\frac{48}{577}\right).$

Note that
$\mu_{1}=\lambda_{1}^{2}~{},~{}\mu_{2}=\lambda_{2}^{2}=\lambda_{1}^{-2}$,
whence also $\mu_{1}\cdot\mu_{2}=1$. Consequently:

$z_{n}^{2}=e^{2}\lambda_{1}^{4n}+2eg\left(-1\right)^{n}\lambda_{1}^{2n}+\left(2ef+g^{2}\right)+2fg\left(-1\right)^{n}\lambda_{1}^{-2n}+f^{2}\lambda_{1}^{-4n}$

Finally, we verify the following five equalities:

  * •

$e^{2}=a^{4}+c^{4}$

  * •

$f^{2}=b^{4}+d^{4}$

  * •

$2eg=4a^{3}b+4c^{3}d$

  * •

$2fg=4ab^{3}+4cd^{3}$

  * •

$2ef+g^{2}=6a^{2}b^{2}+6c^{2}d^{2}-8$

which together imply Equation (1).

## References

  * [1] G. H. Hardy and E. M. Wright.  An Introduction to the Theory of Numbers, pages 191–192.  Oxford, fourth edition, 1975.