Lepton–Nucleus Interactions within Microscopic Approaches

Abstract

This review paper emphasizes the significance of microscopic calculations with quantified theoretical error estimates in studying lepton–nucleus interactions and their implications for electron scattering and accelerator neutrino oscillation measurements. We investigate two approaches: Green’s Function Monte Carlo and the extended factorization scheme, utilizing realistic nuclear target spectral functions. In our study, we include relativistic effects in Green’s Function Monte Carlo and validate the inclusive electron scattering cross section on carbon using available data. We compare the flux-folded cross sections for neutrino-carbon scattering with T2K and MINER$\nu$A experiments, noting the substantial impact of relativistic effects in reducing the theoretical curve strength when compared to MINER$\nu$A data. Additionally, we demonstrate that quantum Monte Carlo-based spectral functions accurately reproduce the quasi-elastic region in electron scattering data and T2K flux-folded cross sections. By comparing results from Green’s Function Monte Carlo and the spectral function approach, which share a similar initial target state description, we quantify errors associated with approximations in the factorization scheme and the relativistic treatment of kinematics in Green’s Function Monte Carlo.

1. Introduction

The electron scattering experimental program at Jefferson laboratory, aimed at investigating nuclear short-range correlations [1,2,3], and the accelerator neutrino program, which will culminate with the completion of DUNE [4], have been the springboard for significant progress in theoretical calculations of lepton–nucleus scattering. Approaches based on empirical effective nucleon–nucleon interactions [5,6,7,8] have been used to study inclusive and semi-inclusive neutrino scattering data in a variety of kinematic setups [9,10,11,12,13,14,15,16,17,18,19,20]. Despite their success, it is still imperative to attain a description of lepton–nucleus scattering from microscopic nuclear dynamics, which assumes that the structure and electroweak properties of atomic nuclei can be modeled in terms of nuclear potentials and consistent electroweak currents. These microscopic approaches allow one to quantify the theoretical uncertainties due to both modeling nuclear dynamics and solving the many-body Schrödinger equation. This aspect is critical for a meaningful comparison with electron scattering data, and, perhaps more importantly, to rigorously assess the error budget of neutrino oscillation parameters. Moreover, retaining nuclear correlations in the initial target state is important to explain the observed abundances of neutron–proton correlated pairs with respect to the proton–proton and neutron–neutron ones [21]. These experimental measurements can in turn shed light on the behavior of nuclear forces at short distances, which plays an important role in the equation of the state of infinite nuclear matter at high density [22,23].

Variational Monte Carlo (VMC) and Green’s Function Monte Carlo (GFMC) methods have proven to be extremely successful for computing the structure and electroweak transitions of atomic nuclei taking as input highly realistic nuclear Hamiltonians [24]. Over the past decade, these methods have been employed to carry out microscopic calculations of the electroweak response functions of light nuclei, fully retaining correlations and consistent one- and two-body currents [25,26,27,28]. Computing the hadronic response tensor is a highly nontrivial task, as it involves transitions from the initial ground state of the target to excited states, both bound and in the continuum. The prohibitive difficulties involved in computing all transitions mediated by the electroweak current operators are circumvented by employing integral transform techniques. Within this approach, the electroweak response functions are inferred from their Laplace transforms, denoted as Euclidean responses, which are estimated during the GFMC imaginary-time propagation. Retrieving the energy dependence of the response functions from their Euclidean counterparts is nontrivial. The maximum entropy method [29,30] has been extensively employed to retrieve the energy dependence of the electroweak response functions in the smooth quasi-elastic region. More recently, inversion methods based on deep neural networks have been proposed as viable alternatives and seem to be more accurate, especially in the low-energy transfer region [31].

One of the main limitations of the GFMC approach lies in the nonrelativistic formulation of the many-body problem. Although the leading relativistic corrections are included in the transition operators [32], the kinematics of the reaction is nonrelativistic, thereby limiting the application of the GFMC to moderate values of the momentum transfer. This restriction is particularly relevant when making predictions for inclusive neutrino–nucleus cross sections since the incoming neutrino flux is not monochromatic and its tails extend to high energies. In Refs. [33,34], relativistic effects in GFMC calculations of lepton–nucleus scattering are controlled by choosing a reference frame that minimizes nucleon momenta and utilizing the so-called “two-fragment” model to include relativity in the kinematics of the reaction.

On the other hand, alternative approaches based on the factorization of the nuclear final state, such as the spectral function (SF) formalism [35], can reach larger energies and momentum transfers, as they include relativistic effects in both the kinematics and in the interaction vertex. In contrast with the GFMC, the SF approach can access exclusive channels and larger nuclei. However, while based on a similar treatment of the initial target state, factorizing the final state involves additional approximations, which are only valid at large momentum transfer and whose validity can be tested against comparisons with GFMC calculations [36,37].

In this work, we first review the GFMC and SF approaches to compute inclusive lepton–nucleus scattering, placing particular emphasis on the role of relativistic effects and two-body currents. We then compare the SF predictions for the neutrino–nucleus cross sections with the MINER$\nu$A Medium-Energy charge-current quasielastic (CCQE)-like data [38]. Finally, we present unpublished GFMC calculations for the inclusive electron-${}^{12}$C cross sections that include relativistic corrections.

2. Methodology

The lepton–nucleus differential cross section in the one-boson exchange approximation can be written as

$${\left(\frac{d\sigma }{d{E}^{\prime }d{\Omega }^{\prime }}\right)}_l={\mathcal{C}}_l{L}_{\mu \nu }^l{R}^{\mu \nu },$$

(1)

where l stands for either a charged lepton or neutrino, $C_l$ is a coupling term, and $E^{\prime }$ and ${\Omega }^{\prime }$ are the energy and solid angle of the lepton in the final state. The leptonic tensor is denoted by $L_{\mu \nu }^l$ and is a function of the initial and final lepton four-momenta k and $k^{\prime }$, respectively. For small lepton energy, the Coulomb distortion of the outgoing lepton in the potential of the residual nucleus can be described by multiplying the cross section by the Fermi function $F(Z,k^{\prime })$, with Z denoting the number of protons. The expression of this function for charge-raising reactions is given in Ref. [32], and it is equal to one otherwise. For higher energies, the correction is provided by the modified effective momentum approximation as discussed in Ref. [39], where an effective momentum is utilized for the final lepton, correcting its value with the Coulomb energy evaluated at the center of the nucleus and modifying the phase space representing the density of final states accordingly. For the comparisons with T2K and MINER$\nu$A data discussed in this review, the effect of Coulomb corrections is negligible as discussed in Figure 9 of Ref. [40] and therefore they have not been included.

In this review, we will consider electromagnetic and CC electroweak interactions. In the first case, we have an electron in both the initial and final state, the prefactor reads $C_{e^{-}}=\alpha E/\left(Q^4{E}^{\prime }\right)$, where E is the energy of the initial lepton, $\alpha =1/137$ is the electromagnetic fine structure constant, $Q^2={\mathbf{q}}^2-{\omega }^2$ is the four-momentum transfer and the leptonic tensor is

$$\begin{array}{c}\hfill L_{\mu \nu }=\frac{1}{E{E}^{\prime }}(k_{\mu }{k}_{\nu }^{\prime }+k_{\mu }^{\prime }{k}_{\nu }-g_{\mu \nu }(k·k^{\prime }-m_e^2)),\end{array}$$

(2)

where $m_e=511$ keV is the electron mass. Note that we adopted the convention $h=c=1$. For CC electroweak interactions, we have that a neutrino or anti-neutrino scatters off the initial nucleus and in the final state the corresponding charge lepton is emitted. The prefactor reads $C_l={(G_{F}cos{\theta }_c)}^2/\left(4\pi \right)\left|{\mathbf{k}}^{\prime }\right|E^{\prime }$ with $G_F=1.1803\times {10}^{-5}$ [41] and $cos{\theta }_c=0.97425$ [42]. The leptonic tensor has an additional term proportional to the Levi-Civita tensor

$$\begin{array}{c}\hfill L_{\mu \nu }=\frac{1}{E{E}^{\prime }}(k_{\mu }{k}_{\nu }^{\prime }+k_{\mu }^{\prime }{k}_{\nu }-g_{\mu \nu }k·k^{\prime }\pm {ϵ}_{\mu \rho \nu \sigma }{k}^{\rho }{k^{\prime }}^{\sigma }),\end{array}$$

(3)

where the sign +(-) corresponds to a $\nu$ ($\overline{\nu }$) in the initial state.

The hadronic response tensor, $R^{\mu \nu }$, contains all the information on the structure of the nuclear target and is defined as

$$R_{\mu \nu }=\sum _f\langle {\mathrm{\Psi }}_0|J_{\mu }^{†}|{\mathrm{\Psi }}_f\rangle \langle {\mathrm{\Psi }}_f|J_{\nu }|{\mathrm{\Psi }}_0\rangle \delta (E_0+\omega -E_f),$$

(4)

in terms of a sum over all transitions from the ground state $|{\mathrm{\Psi }}_0\rangle$ with energy $E_0$ to any final state $|{\mathrm{\Psi }}_f\rangle$ with energy $E_f$, including states with additional hadrons. The nuclear current operator describing the interaction with the electroweak probe is denoted by $J_{\mu }$.

2.1. Nuclear Hamiltonian and Current Operator

Microscopic nuclear methods are aimed at describing properties of nuclear systems as they emerge from the individual interactions among the constituent protons and neutrons. This endeavor is based on the tenet that the internal structure of atomic nuclei can be described starting from a nonrelativistic Hamiltonian of A point-like nucleons

$$H=\sum _i^A\frac{{\mathbf{p}}_i^2}{2m_N}+\sum _{i<j}^A{v}_{ij}+\sum _{i<j<k}^A{V}_{ijk}.$$

(5)

In the above equation, $\mathbf{p}$ and $m_N$ are the nucleon momentum and mass defined as the average of the proton and neutron mass $m_N=(m_p+m_n)/2$, while $v_{ij}$ and $V_{ijk}$ are the two- (NN) and three-nucleon (3N) potentials, respectively; four- and higher-body potentials are assumed to be suppressed.

Phenomenological NN interactions have been traditionally constructed by including the long-range one-pion exchange interaction, while different schemes are implemented to account for intermediate- and short-range effects, including multiple-pion exchange, contact terms, heavy-meson exchange, or the excitation of nucleons into virtual $\mathsf{\Delta }$-isobars. As an example, the highly accurate Argonne $v_{18}$ (AV18) potential [43] involves a number of parameters that are determined by fitting deuteron properties and the large database of NN scattering data at laboratory energies up to pion production thresholds. The AV18 potential is written as

$$v_{ij}=\sum _{p=1}^{18}{v}_p\left(r_{ij}\right)O_{ij}^p.$$

(6)

The first 14 spin–isospin operators are charge-independent

$$O_{ij}^{1-14}=[1,{\sigma }_i·{\sigma }_j,S_{ij},\mathbf{L}·\mathbf{S},{\mathbf{L}}^2,{\mathbf{L}}^2({\sigma }_i·{\sigma }_j),{(\mathbf{L}·\mathbf{S})}^2]]\times [1,{\tau }_i·{\tau }_j],$$

(7)

where ${\sigma }_i$ are Pauli matrices that operate over the spin of nucleons, $S_{ij}=3({\widehat{r}}_{ij}·{\sigma }_i)({\widehat{r}}_{ij}·{\sigma }_j)-{\sigma }_i·{\sigma }_j$ is the tensor operator, ${\mathbf{L}}_{ij}=\frac{1}{2i}({\mathbf{r}}_i-{\mathbf{r}}_j)\times ({\nabla }_i-{\nabla }_j)$ is the relative angular momentum of the pair $ij$ and $\mathbf{S}=\frac{1}{2}({\nabla }_i+{\nabla }_j)$ is the total spin. The remaining operators include three charge-dependent terms and one charge-symmetry-breaking contribution

$$\begin{array}{cc}\hfill O_{ij}^{15-17}& =[1,{\mathbf{\sigma }}_i·{\mathbf{\sigma }}_j,S_{ij}]\times T_{ij},\hfill \\ \hfill O_{ij}^{18}& ={\tau }_i^z+{\tau }_j^z,\hfill \end{array}$$

(8)

where $T_{ij}=3{\tau }_i^z{\tau }_j^z-{\mathbf{\tau }}_i·{\mathbf{\tau }}_j$ is the isotensor operator.

Phenomenological 3N interactions, consistent with the NN ones, are generally expressed as a sum of a two-pion-exchange P-wave term, a two-pion-exchange S-wave contribution, a three-pion-exchange contribution, plus a contact interaction. Their inclusion is essential for reproducing the energy spectrum of atomic nuclei and the saturation properties of infinite nucleonic matter. For instance, the Illinois-7 3N force [44], when used together with AV18, can reproduce the spectrum of nuclei up to ${}^{12}$C with percent-level accuracy—see Figure 2 discussed in Section 3.

The past two decades have witnessed the tremendous development and success of chiral Effective Field Theory [45,46,47,48,49,50,51,52,53,54,55,56] ($\chi$EFT). This formalism exploits the broken chiral symmetry pattern of QCD, the fundamental theory of strong interactions, to construct an effective Hamiltonian organized in powers of the ratio between the pion mass, $m_{\pi }$, or a typical nucleon momentum, Q, and the scale of chiral symmetry breaking, ${\mathrm{\Lambda }}_{\chi }\sim 1$ GeV. Over the years, NN interactions have been developed up to ${\mathrm{N}}^5\mathrm{LO}$ in the chiral expansion [57,58,59], with a full systematic error analysis currently underway [60]. On the other hand, chiral 3N forces have been fully derived at ${\mathrm{N}}^3\mathrm{LO}$, while only contact terms at ${\mathrm{N}}^4\mathrm{LO}$ have so far been included [61].

Analogously to the nuclear Hamiltonian, the nuclear current operator $J^{\mu }$, which couples the nucleus to the external electroweak probe, can be written as a sum of both one- and two-body contributions

$$J_A^{\mu }\left(q\right)=\sum _i{j}_i^{\mu }\left(q\right)+\sum _{ij}{j}_{ij}^{\mu }\left(q\right)+\dots$$

(9)

where higher-order terms, involving three nucleons or more, are found to be small [62] and generally neglected.

The one-body electromagnetic current is given by

$$\begin{array}{c}\hfill j_{\mathrm{EM}}^{\mu }\left(q\right)=j_{\gamma ,S}^{\mu }\left(q\right)+j_{\gamma ,z}^{\mu }\left(q\right),\end{array}$$

(10)

where the first term is the isoscalar contribution and the second one is the isovector. The isoscalar component reads

$$\begin{array}{c}\hfill j_{\gamma ,S}^{\mu }\left(q\right)=\frac{G_E^S+\tau G_M^S}{2(1+Q^2/4m_N^2)}{\gamma }^{\mu }+i\frac{{\sigma }^{\mu \nu }{q}_{\nu }}{4m_N}\frac{G_M^S-G_E^S}{1+Q^2/4m_N^2}.\end{array}$$

(11)

The isoscalar and isovector components of electric and magnetic form factors are written in terms of the proton and neutron ones as

$$\begin{array}{c}\hfill G_{E,M}^S=G_{E,M}^p+G_{E,M}^n,G_{E,M}^V=G_{E,M}^p-G_{E,M}^n.\end{array}$$

(12)

The isovector contribution to the current operator $j_{\gamma ,z}^{\mu }$ is obtained by replacing $G_{E,M}^S\to G_{E,M}^V{\tau }_z$ in Equation (11).

The one-body charge and current operators employed in the GFMC are obtained from the nonrelativistic reduction of the covariant operator of Equation (11), including all the terms up to $1/m_N^2$. This expansion leads to the following expressions for isoscalar charge, both transverse (⊥) and longitudinal (‖) to $\mathbf{q}$ components of the current operator

$$\begin{array}{cc}\hfill j_{\gamma ,S}^0=& \frac{G_E^S}{2\sqrt{1+Q^2/4m_N^2}}-i\frac{2G_M^S-G_E^S}{8m_N^2}\mathbf{q}·(\mathbf{\sigma }\times \mathbf{p})\hfill \\ \hfill {\mathbf{j}}_{\gamma ,S}^{\perp }=& \frac{G_E^S}{2m_N}{\mathbf{p}}^{\perp }-i\frac{G_M^S}{4m_N}(\mathbf{q}\times \mathbf{\sigma })\hfill \\ \hfill j_{\gamma ,S}^{\Vert }=& \frac{\omega }{\left|\mathbf{q}\right|}{j}_{\gamma ,S}^0.\hfill \end{array}$$

(13)

Note that the last relation has been obtained from the conserved vector current (CVC) relation [63], e.g., $\omega J^0(\omega ,\mathbf{q})-\mathbf{q}·\mathbf{J}(\omega ,\mathbf{q})=0$. The CC electroweak interactions of a neutrino or anti-neutrino with the hadronic target are written as the sum of a vector and axial term

$$\begin{array}{cc}\hfill j_{CC}^{\mu }\left(q\right)=& j_{V,\pm }^{\mu }\left(q\right)+j_{A,\pm }^{\mu }\left(q\right).\hfill \end{array}$$

(14)

The CVC hypothesis allows one to write $j_{V,\pm }^{\mu }\left(q\right)$ in terms of the isovector term where ${\tau }_z$ is replaced by the isospin raising–lowering operator ${\tau }_{\pm }=({\tau }_{i,x}\pm {\tau }_{i,y})/2$. The relativistic expression of the axial one-body current operator reads

$$\begin{array}{c}\hfill j_{A,\pm }^{\mu }=-{\gamma }^{\mu }{\gamma }_5{G}_A{\tau }_{\pm }-q^{\mu }{\gamma }_5\frac{G_P}{m_N}{\tau }_{\pm }.\end{array}$$

(15)

Based on Partially Conserved Axial Current (PCAC) arguments, the pseudo-scalar form factor is written in terms of the axial one

$$\begin{array}{c}\hfill G_P=\frac{2m_N^2}{(m_{\pi }^2+Q^2)}{G}_A.\end{array}$$

(16)

Most neutrino–nucleus scattering calculations are carried out employing a dipole parameterization for the axial form factor, which is given by

$$\begin{array}{cc}\hfill G_A& =\frac{g_A}{{(1+Q^2/{\mathrm{\Lambda }}_A^2)}^2},\hfill \end{array}$$

(17)

where the nucleon axial-vector coupling constant is taken to be $g_A=1.2694$ [64] and the axial mass is taken as ${\mathrm{\Lambda }}_A=1.049$ GeV [65]. More recently, a model-independent z expansion has been introduced to parameterize the axial form factor

$$G_A\left(Q^2\right)=\sum _{j=0}^{\infty }{a}_{j}z{\left(Q^2\right)}^j\approx \sum _{j=0}^{j_{\max }}{a}_{j}z{\left(Q^2\right)}^j.$$

(18)

In the last equation, z is an analytic function of $Q^2$ for $Q^2=-t>-t_c$

$$z\left(Q^2\right)=\frac{\sqrt{t_c+Q^2}-\sqrt{t_c-t_0}}{\sqrt{t_c+Q^2}+\sqrt{t_c-t_0}},$$

(19)

where $t_c$ is the location of the t-channel cut [66,67,68] and $t_0$ is an arbitrary parameter. The coefficient $a_j$ includes nucleon structure information and $j_{\max }$ is a truncation parameter required to make the number of expansion parameters finite. The coefficients of this expansion are determined by fitting either neutrino–deuteron scattering data [69] or Lattice-QCD nucleon axial current matrix elements at several discrete values of $Q^2$ [70,71,72]. The results of these different determinations of the axial form factors are displayed in Figure 1. While an agreement between different LQCD calculations is clearly visible, the LQCD axial form factor results are 2–3$\sigma$ larger than the results of Ref. [69] for $Q^2\gtrsim 0.3{\mathrm{GeV}}^2$. The impact of these tensions in the $Q^2$ dependence of the axial form factor on neutrino–nucleus cross section predictions has been discussed in Refs. [73,74] and recently in [37,75].

For the CC processes, we report the nonrelativistic reduction of the charge and axial current operators [77] (for brevity, we neglect order $1/m_N^2$ terms)

$$\begin{array}{c}\hfill j_{A,\pm }^0=-\frac{G_A}{2m_N}{\tau }_{\pm }\mathbf{\sigma }·(\mathbf{q}+\mathbf{p}),{\mathbf{j}}_{A,\pm }=-\frac{G_A}{2m_N}\mathbf{\sigma }{\tau }_{\pm }.\end{array}$$

(20)

and the pseudoscalar contribution

$$\begin{array}{c}\hfill j_{PS,\pm }^{\mu }=\frac{G_A}{m_{\pi }^2+Q^2}{\tau }_{\pm }{q}^{\mu }\mathbf{\sigma }·\mathbf{q}.\end{array}$$

(21)

The current conservation relation can be rewritten as

$$\nabla ·{\mathbf{J}}_{\mathrm{EM}}+i[H,J_{\mathrm{EM}}^0]=0.$$

(22)

It requires the introduction of a two-body current operator in $J_{\mathrm{EM}}^{\mu }$ and links the divergence of this operator to the commutator of the charge operator with the nucleon–nucleon interaction.

For the electromagnetic case, however, gauge invariance actually puts constraints on these form factors by linking the divergence of the two-body currents to the commutator of the charge operator with the nucleon–nucleon interaction.

Within $\chi$EFT, one can exploit the gauge invariance of the theory and construct nuclear current operators that are fully consistent with the nuclear potentials, at each order of the chiral expansion. The derivation of $\chi$EFT two-body electroweak currents has been the subject of extensive studies carried out by different groups [78,79,80,81,82,83,84].

The majority of the results that will be presented in this review have been obtained utilizing semi-phenomenological currents that are consistent with the AV18 potential. The isoscalar and isovector components of the two-body electromagnetic current operator consist of “model-independent” and “model-dependent term” terms. The former are obtained from the NN interaction, and by construction satisfy current conservation. They consist of the one-pion and one-rho exchange current operators—their expressions are well known and reported in Refs. [35,85,86], both in their relativistic and non relativistic formulation.

The transverse components of the two-body currents cannot be directly linked to the nuclear Hamiltonian. The isovector current is associated with the exchange of a pion followed by the excitation of a $\mathsf{\Delta }$-resonance in the intermediate state. The isoscalar contribution includes the $\rho \pi \gamma$ transition, whose couplings are extracted from the widths of the radiative decay $\rho \to \pi \gamma$ and the $Q^2$ dependence of the electromagnetic transition form factor is modeled assuming vector-meson dominance [77,87].

3. Quantum Monte Carlo Approaches

Solving the Schrödinger equation for the nuclear Hamiltonian defined in Equation (5) entails nontrivial difficulties, owing to the nonperturbative nature and strong spin–isospin dependence of realistic nuclear forces. The VMC method is routinely employed to approximately find the ground-state solution of the quantum many-body problem for nuclei with up to $A=12$ nucleons [24]. Within this approach, the true ground state ${\mathrm{\Psi }}_0$ is approximated by a variational state ${\mathrm{\Psi }}_V$, which is defined in terms of a set of variational parameters. The optimal values of the latter are found by exploiting the variational principle, i.e., by minimizing the variational energy

$$E_V=\frac{\langle {\mathrm{\Psi }}_V\left|H\right|{\mathrm{\Psi }}_V\rangle }{\langle {\mathrm{\Psi }}_V|{\mathrm{\Psi }}_V\rangle }\ge E_0.$$

(23)

The form of the variational state is taken to be

$${\mathrm{\Psi }}_V=\mathcal{F}|\mathrm{\Phi }\rangle ,$$

(24)

where $\mathcal{F}$ is a permutation-invariant correlation operator of a Jastrow, and the anti-symmetric $|\mathrm{\Phi }\rangle$ controls the quantum numbers and the long-range behavior of the wave function. The correlation operator explicitly includes correlations between pairs and triplets of nucleons

$$\mathcal{F}=\left(\mathcal{S}\prod _{i<j<k}(1+F_{ijk})\right)\left(\mathcal{S}\prod _{i<j}{F}_{ij}\right),$$

(25)

where S is the symmetrization operator, which is required to ensure the anti-symmetry of ${\mathrm{\Psi }}_V$ since, in general, neither the two-body correlations, $F_{ij}$, nor the three-body ones, $F_{ijk}$, commute. The structure of the spin-dependent nuclear correlation operators reflects the one of the NN potential of Equation (6)

$$F_{ij}=\sum _{p=1}^6{u}^p\left(r_{ij}\right)O_{ij}^p,$$

(26)

where the first six operators of Equation (7) are $O_{ij}^{1-6}=[1,{\mathbf{\sigma }}_i·{\mathbf{\sigma }}_j,S_{ij},]\times [1,{\mathbf{\tau }}_i·{\mathbf{\tau }}_j]$. More sophisticated correlation operators that explicitly include spin–orbit correlations have been used in the cluster-variational Monte Carlo calculations of Ref. [88]. However, the computational cost of these additional terms is significant, while the the gain in the variational energy is relatively small [89].

The GFMC evolves the variational state in imaginary time to filter out the excited state components so that

$$|{\mathrm{\Psi }}_0\rangle =\underset{\tau \to \infty }{lim}|\mathrm{\Psi }\left(\tau \right)\rangle =\underset{\tau \to \infty }{lim}\exp [-(H-E_0)\tau ]|{\mathrm{\Psi }}_T\rangle ,$$

(27)

The above imaginary-time evolution is carried out as a series of many small steps $\mathsf{\Delta }\tau$ using an exact two-body short-time propagator [90]. At each step, the GFMC retains all of the spin–isospin components of the nuclear wave function and can take as input the most realistic local interactions. The results for the ground-state energies of nuclei up to ${}^{12}$C has been computed with 1% accuracy within GFMC using the semi-phenomenological AV18+IL7 potentials in Ref. [24], and they are displayed in Figure 2. Note that a plot with a comparable degree of accuracy has also been obtained, using as input the $\mathsf{\Delta }$-full $\chi$EFT nuclear forces that are local in coordinate space [91,92].

Figure 2. Energies of light nuclear ground and excited states from a particular parameterization of Equation (5) computed using Green’s Function Monte Carlo (GFMC) techniques. Figure from [24].

Figure 2. Energies of light nuclear ground and excited states from a particular parameterization of Equation (5) computed using Green’s Function Monte Carlo (GFMC) techniques. Figure from [24].

Since all the spin–isospin degrees of freedom are retained, the GFMC suffers from an exponential scaling with the number of nucleons, which currently limits its applicability to light nuclei, up to ${}^{12}$C. The Auxiliary Field Diffusion Monte Carlo can reach larger nuclear systems by representing the spin–isospin degrees of freedom in terms of products of single-particle states, thereby reducing the computational cost from exponential to polynomial in A [24,93]. However, the use of Hubbard–Stratonovich transformations in the AFDMC imaginary-time propagation prevents the AFDMC from treating highly realistic NN potentials that include an isospin-dependent spin–orbit term.

3.1. Green’s Function Monte Carlo Calculations of Electroweak Responses

GFMC techniques go beyond just the calculation of ground-state energies and wave functions. The dynamical properties of the nucleus can be extracted by reducing the sum over the final states in Equation (4) to the expectation value of a kernel operator evaluated in the ground state. More specifically, we consider the Euclidean response function

$$\begin{array}{cc}\hfill E_{\alpha \beta }(\mathbf{q},\tau )& =\int d\omega K(\tau ,\omega )R_{\alpha \beta }(\mathbf{q},\omega )\hfill \\ \hfill & =\sum _f\langle {\mathrm{\Psi }}_0|J_{\alpha }^{†}\left(\mathbf{q}\right)|{\mathrm{\Psi }}_f\rangle K(\tau ,E_f-E_0)\langle {\mathrm{\Psi }}_f|J_{\beta }\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle ,\hfill \end{array}$$

where $K(\tau ,\omega )$ is a yet-to-be-specified kernel. Using a completeness relation amongst the final states, this can be simplified to

$$E_{\alpha \beta }(\mathbf{q},\tau )=\sum _f\langle {\mathrm{\Psi }}_0|J_{\alpha }^{†}\left(\mathbf{q}\right)K(\tau ,H-E_0)J_{\beta }\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle$$

(28)

so that the problem involves only the ground state. Choosing an appropriate kernel function allows one to solve for the Euclidean response using ab initio methods. In particular, a Laplace kernel has been adopted with GFMC techniques, yielding the following expression for the inelastic contribution to the response function

$$\begin{array}{cc}\hfill E_{\alpha \beta }(\mathbf{q},\tau )& ={\int }_0^{\infty }d\omega R_{\alpha \beta }{e}^{-\omega \tau }=\langle {\mathrm{\Psi }}_0|J_{\alpha }^{†}\left(\mathbf{q}\right)e^{-(H-E_0)\tau }{J}_{\beta }\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle .\hfill \end{array}$$

(29)

In the electromagnetic case, only the longitudinal $(R_L=R_{00})$ and transverse $(R_T=R_{xx}+R_{yy})$ responses contribute. In the longitudinal case, we remove the elastic contribution, in which the final state is simply the recoiling ground state, by defining

$$\begin{array}{cc}\hfill E_{00}(\mathbf{q},\tau )& ={\int }_{{\omega }_{el}^{+}}^{\infty }d\omega R_{00}{e}^{-\omega \tau }=\langle {\mathrm{\Psi }}_0|J_0^{†}\left(\mathbf{q}\right)e^{-(H-E_0)\tau }{J}_0\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle -|F_0{\left(\mathbf{q}\right)|}^2{e}^{-{\omega }_{el}\tau },\hfill \end{array}$$

(30)

In the above equation, ${\omega }_{el}={\mathbf{q}}^2/2M_A$, with $M_A$ being the mass of the nucleus, is the energy of the recoiling ground state and the elastic form factor is defined as $F_0\left(\mathbf{q}\right)=\langle {\mathrm{\Psi }}_0|J_0\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle$.

The calculation of the imaginary-time correlation operator in the right hand side of Equation (30) follows the same methodology applied to project out the exact ground state of H from a trial wave function in Equation (27). First, an unconstrained imaginary-time propagation of the state $|{\mathrm{\Psi }}_0\rangle$ is performed and stored. Then, the states $J_{\alpha }\left(\mathbf{q}\right)|{\mathrm{\Psi }}_0\rangle$ are evolved in imaginary time following the path previously saved. For a complete discussion of the methods, see Refs. [27,94,95]. To retrieve the energy dependence of the response functions, Bayesian techniques, most notably maximum entropy (MaxEnt), have been developed specifically for this type of problem [94] and successfully exploited to obtain smooth quasi-elastic responses [27,95]. However, MaxEnt struggles to reconstruct the narrow peaks corresponding to low-energy transitions. In particular, understanding the low-lying nuclear transitions is necessary to properly describe the longitudinal electromagnetic responses of ${}^{12}$C in low-energy transfer. The results of Ref. [27] have been obtained by subtracting the contribution of these excited states by defining

$$\begin{array}{c}\hfill {\overline{E}}_{00}(\mathbf{q},\tau )=E_{00}(\mathbf{q},\tau )-\sum _f{\left|\langle {\mathrm{\Psi }}_f\left|J_0\left(\mathbf{q}\right)\right|{\mathrm{\Psi }}_0\rangle \right|}^2{e}^{-(E_f-E_0)\tau },\end{array}$$

(31)

where the sum only includes the $2^{+}$, $0_2^{+}$ and $4^{+}$ final states. The experimental energies and longitudinal transition form factors from Refs. [96,97] are used.

Furthermore, understanding this region is also crucial for detecting supernova neutrinos as well as to describe the low-energy tail contribution of the neutrino flux in accelerator experiments. To this aim, in Ref. [31] an exploratory study was carried out to develop physics-informed artificial neural network architectures suitable for approximating the inverse of the Laplace transform, utilizing simulated, albeit realistic, electromagnetic response functions. The training was performed using pairs of physically meaningful responses and their Laplace transform. There are two data sets of response functions, characterized by either one or two distinct peaks in the energy-transfer domain. The left panel of Figure 3 displays a subset of the two-peaks training data. A detailed comparison of the reconstruction results obtained for both the one- and two-peak data sets demonstrates that the physics-informed artificial neural network outperforms MaxEnt in both the low-energy transfer and the quasi-elastic regions—an illustrative example of this trend is shown in the right panel of Figure 3. Work is currently underway to extend the study of Ref. [31] to real GFMC data and to perform error propagation.

3.1.1. Relativistic Corrections

One of the limitations of the GFMC approach to describe nuclear reactions is the nonrelativistic formulation of the many-body problem. Although the leading relativistic corrections are typically included in the transition operators [32], the kinematics of the reaction is treated as nonrelativistic, and an expansion of fully relativistic currents in $p/m$ is performed. The explicit expressions of the one-body current operators adopted in the GFMC calculation are reported in Section 2.1. Thereby, the application of these methods is limited to moderate values of the momentum transfer.

In a number of studies [33,34,98,99,100,101,102], a method has been proposed to extend the applicability of manifestly nonrelativistic hyperspherical harmonics and Quantum Monte Carlo (QMC) methods to higher momentum-transfer values than typically possible. This method reduces relativistic effects by performing the calculations in a reference frame that minimizes nucleon momenta. The reference frame that achieves this goal for kinematics close to the quasi-elastic peak is the active nucleon Breit frame (ANB). The ANB is defined as the reference frame moving along the direction of the momentum transfer $\mathbf{q}$ where ${\mathbf{P}}_i^{ANB}=-A{\mathbf{q}}^{ANB}/2$, with ${\mathbf{P}}_i^{ANB}$ the momentum of the initial nucleus in the ANB. Indeed, if one assumes that the bulk of the momentum is transferred to a single nucleon, in the ANB this nucleon has initial momentum $\mathbf{k}\approx {\mathbf{P}}_0^{ANB}/A=-{\mathbf{q}}^{ANB}/2$. The corresponding final-state nucleon has momentum $\mathbf{k}+{\mathbf{q}}^{ANB}={\mathbf{q}}^{ANB}/2$. Hence, in the ANB the magnitude of both the initial and final-state nucleon momentum is minimal. Additionally, the energy transfer at the quasi-elastic peak is zero in the ANB frame, implying that ${\mathbf{q}}^{ANB}$ is also minimal at the quasi-elastic peak compared to other frames.

Within a nonrelativistic calculation, the nuclear response can be computed in different reference frames by evaluating Equation (4) at the momentum transfer in the reference frame specified by ${\mathbf{q}}^{fr}$, and by taking into account the kinetic energy of initial- and final-state systems in the energy balance. Thus, the energy-conserving delta function $\delta (E_0-\omega -E_f)$ is evaluated with $\omega ={\omega }^{fr}$, and $E_0-E_f=-{\left({\mathbf{P}}_f^{fr}\right)}^2/\left(2M_A\right)+{\left({\mathbf{P}}_0^{fr}\right)}^2/\left(2M_A\right)$, leading to an energy shift of the response.

The dependence on the reference frame used for calculations can be evaluated by performing a Lorentz boost of the response back to the LAB frame. At momentum transfers larger than 500 $\mathrm{MeV}$, one starts to see differences between calculations performed in different reference frames [33,34,98], indicating that relativistic effects become important.

This frame dependence in the region of the quasi-elastic peak can be significantly reduced by including the assumption of single-nucleon knockout in the energy balance. In order to achieve this, one can use the so-called two-fragment model, where a breakup into two fragments, the nucleon and residual system, is assumed. Following the arguments of Refs. [33,98], the approach consists of evaluating the nuclear response at an energy $p_{rel}^2/\left(2\mu \right)$, with $p_{rel}$ being the magnitude of the relative momentum of the two fragments and $\mu$ being the reduced mass. The energy of the final-state system can be written in a relativistic way as

$$\begin{array}{cc}\hfill \omega +E_0=E_f& =\sqrt{m^2+{({\mathbf{p}}_{rel}+(\mu /M){\mathbf{P}}_f)}^2}\hfill \\ \hfill & +\sqrt{M^2+{({\mathbf{p}}_{rel}-(\mu /m){\mathbf{P}}_f)}^2};\hfill \end{array}$$

(32)

where ${\mathbf{P}}_f={\mathbf{P}}_0+\mathbf{q}$ is the center-of-mass momentum. Under the assumption that ${\mathbf{p}}_{rel}$ is directed along $\mathbf{q}$, one can solve Equation (32) for $p_{rel}$.

In Refs. [33,98], it is indeed found that the frame dependence for electroweak scattering is strongly reduced when including the two-fragment model to determine the energy. Moreover, the resulting LAB frame responses are practically identical to the response obtained in the ANB when the fragment model is not included [34].

Calculations of the nuclear response in the ANB can be used to extend the applicability of GFMC responses to larger momentum transfer. In Ref. [33], an improved description of $(e,e^{\prime })$ data for scattering off ${}^4$He was obtained at large momentum transfer with GFMC responses computed in the ANB. Recently, this approach was applied to GFMC calculations of flux-folded charged-current neutrino scattering off ${}^{12}$C [34].

4. Extended Factorization Scheme

At large values of the momentum transfer, ($\left|\mathbf{q}\right|\gtrsim 400$ MeV), the Impulse Approximation (IA) can be applied in which the lepton–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.

4.1. One-Body Currents

We begin with retaining only one-body current terms and factorize the final state according to

$$|{\mathrm{\Psi }}_f\rangle =|{\mathbf{p}}^{\prime }\rangle \otimes |{\mathrm{\Psi }}_f^{A-1},{\mathbf{p}}_{A-1}\rangle ,$$

(33)

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

$$\langle {\mathrm{\Psi }}_f|j^{\mu }|{\mathrm{\Psi }}_0\rangle \to \sum _k[\langle {\mathrm{\Psi }}_f^{A-1}|\otimes \langle k\left|\right]|{\mathrm{\Psi }}_0\rangle \langle p|\sum _i{j}_i^{\mu }|k\rangle ,$$

(34)

where $\mathbf{p}=\mathbf{q}+\mathbf{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_i^{\mu }$ are specified as the single-nucleon states are just free Dirac spinors. It is important to point out that factorization allows for an account of relativistic effect by adopting Dirac quadri-spinors for the description of the struck particles in the initial and final states and the current operator of Equation (11). These effects become extremely important at large values of $\mathbf{q}$ and $\omega$ where a nonrelativistic calculation is no longer reliable. Substituting the last equation into Equation (4), and exploiting momentum conservation at the single nucleon vertex, allows us to rewrite the incoherent contribution to the one-body hadron tensor as

$$\begin{array}{cc}\hfill R_{1b}^{\mu \nu }(\mathbf{q},\omega )=& \int \frac{d^{3}k}{{\left(2\pi \right)}^3}dE{P}_h(\mathbf{k},E)\frac{m_N^2}{e\left(\mathbf{k}\right)e(\mathbf{k}+\mathbf{q})}\times \sum _i\langle k|j_i^{\mu †}|k+q\rangle \langle k+q|j_i^{\nu }|k\rangle \hfill \\ \hfill & \times \delta (\tilde{\omega }+e\left(\mathbf{k}\right)-e\left(\mathbf{p}\right)),\hfill \end{array}$$

(35)

where $e\left(\mathbf{k}\right)=\sqrt{m_N^2+{\mathbf{k}}^2}$. The factors $m_N/e\left(\mathbf{k}\right)$ and $m_N/e(\mathbf{k}+\mathbf{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\left(\mathbf{k}\right)$ to account for scattering off of a bound nucleon. Finally, the calculation of the one-nucleon spectral function $P_h(\mathbf{k},E)$ provides the probability of removing a nucleon with momentum $\mathbf{k}$ and leaving the residual nucleus with an excitation energy E; its derivation will be discussed in Section 5.

4.2. Two-Body Currents

To describe amplitudes including two-nucleon currents, the factorization ansatz of Equation (33) can be generalized as

$$\begin{array}{c}\hfill |{\mathrm{\Psi }}_f\rangle \to |p{p}^{\prime }{\rangle }_a\otimes |{\mathrm{\Psi }}_f^{A-2}\rangle .\end{array}$$

(36)

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. [35,103,104], the pure two-body current component of the response tensor can be written as

$$\begin{array}{cc}\hfill & R_{2b}^{\mu \nu }(\mathbf{q},\omega )=\frac{V}{2}\int dE\frac{d^{3}k}{{\left(2\pi \right)}^3}\frac{d^3{k}^{\prime }}{{\left(2\pi \right)}^3}\frac{d^{3}p}{{\left(2\pi \right)}^3}\frac{m_N^4}{e\left(\mathbf{k}\right)e\left({\mathbf{k}}^{\prime }\right)e\left(\mathbf{p}\right)e\left({\mathbf{p}}^{\prime }\right)}\hfill \\ \hfill & \times P_h(\mathbf{k},{\mathbf{k}}^{\prime },E)\sum _{ij}\langle k{k}^{\prime }|{j_{ij}^{\mu }}^{†}|p{p}^{\prime }{\rangle }_a\langle p{p}^{\prime }|j_{ij}^{\nu }|k{k}^{\prime }\rangle \delta (\omega -E+2m_N-e\left(\mathbf{p}\right)-e\left({\mathbf{p}}^{\prime }\right)).\hfill \end{array}$$

(37)

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 ${}^{12}$C is taken to be $k_F=225$ MeV. In previous calculations of the above two-body hadron tensor, the two-nucleon spectral function $P_h(\mathbf{k},{\mathbf{k}}^{\prime },E)$ has been approximated as a product of two one-nucleon spectral functions (see Section 5 for a more detailed discussion). This is correct in the infinite nuclear matter limit where the two-nucleon momentum distribution can be split according to

$$n(k,k^{\prime })=n\left(k\right)n\left(k^{\prime }\right)+\mathcal{O}(1/A).$$

(38)

Going beyond this approximation for medium-mass nuclei involves the full calculation of the two-nucleon spectral function including all correlations and will be discussed further in Section 5. The two-body current operator in Equation (37) is given by a sum of four distinct contributions, namely the pion in flight, seagull, pion-pole, and $\mathsf{\Delta }$ excitations

$$\begin{array}{c}\hfill j^{\mu }=\left(j_{\mathrm{pif}}^{\mu }\right)+\left(j_{\mathrm{sea}}^{\mu }\right)+\left(j_{\mathrm{pole}}^{\mu }\right)+\left(j_{\mathsf{\Delta }}^{\mu }\right)\end{array}$$

(39)

and dubbed as Meson-Exchange Currents (MECs). Detailed expressions for each term in Equation (39) can be found in Refs. [35,105]. Below, we only report the two-body current terms involving a $\mathsf{\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 Section 2.1. The current operator can be written as follows [105,106]:

$$\begin{array}{cc}\hfill {(j_{\mathsf{\Delta }}^{\mu })}_{\mathrm{C}C}& =\frac{3}{2}\frac{f_{\pi NN}{f}^{*}}{m_{\pi }^2}\left\{\right[\left(-\frac{2}{3}{\tau }_{\pm }^{\left(2\right)}+\frac{{({\tau }^{\left(1\right)}\times {\tau }^{\left(2\right)})}_{\pm }}{3}\right)F_{\pi NN}\left(k_{\pi }^{\prime }\right)F_{\pi N\mathsf{\Delta }}\left(k_{\pi }^{\prime }\right){\left(j_a^{\mu }\right)}_{\left(1\right)}\hfill \\ \hfill & -\left(\frac{2}{3}{\tau }_{\pm }^{\left(2\right)}+\frac{{({\tau }^{\left(1\right)}\times {\tau }^{\left(2\right)})}_{\pm }}{3}\right)F_{\pi NN}\left(k_{\pi }^{\prime }\right)F_{\pi N\mathsf{\Delta }}\left(k_{\pi }^{\prime }\right){\left(j_b^{\mu }\right)}_{\left(1\right)}]\mathrm{\Pi }{\left(k_{\pi }^{\prime }\right)}_{\left(2\right)}+(1\leftrightarrow 2)\},\hfill \end{array}$$

(40)

where $k^{\prime }$ and $p^{\prime }$ are the initial and final momentum of the second nucleon, respectively, while $k_{\pi }^{\prime }=p^{\prime }-k^{\prime }$ is the momentum of the $\pi$ exchanged in the two depicted diagrams of Figure 4, $f^{*}$ = 2.14, and

$$\begin{array}{cc}\hfill & \mathrm{\Pi }\left(k_{\pi }\right)=\frac{{\gamma }_5{\cancel{k}}_{\pi }}{k_{\pi }^2-m_{\pi }^2},\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & F_{\pi N\mathsf{\Delta }}\left(k_{\pi }\right)=\frac{{\mathrm{\Lambda }}_{\pi N\mathsf{\Delta }}^2}{{\mathrm{\Lambda }}_{\pi N\mathsf{\Delta }}^2-k_{\pi }^2},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & F_{\pi NN}\left(k_{\pi }\right)=\frac{{\mathrm{\Lambda }}_{\pi }^2-m_{\pi }^2}{{\mathrm{\Lambda }}_{\pi }^2-k_{\pi }^2},\hfill \end{array}$$

(43)

with ${\mathrm{\Lambda }}_{\pi N\mathsf{\Delta }}=1150$ MeV and ${\mathrm{\Lambda }}_{\pi }=1300$ MeV. In Equation (40), $j_a^{\mu }$ and $j_b^{\mu }$ denote the $N\to \mathsf{\Delta }$ transition vertices of diagram (a) and (b) of Figure 4, respectively. The expression of $j_a^{\mu }$ is given by

$$\begin{array}{cc}\hfill j_a^{\mu }& ={\left(j_a^{\mu }\right)}_V+{\left(j_a^{\mu }\right)}_A,\hfill \\ \hfill {\left(j_a^{\mu }\right)}_V& ={\left(k_{\pi }^{\prime }\right)}^{\alpha }{G}_{\alpha \beta }\left(p_{\mathsf{\Delta }}\right)[\frac{C_3^V}{m_N}\left(g^{\beta \mu }\cancel{q}-q^{\beta }{\gamma }^{\mu }\right)+\frac{C_4^V}{m_N^2}\left(g^{\beta \mu }q·p_{\mathsf{\Delta }}-q^{\beta }{p}_{\mathsf{\Delta }}^{\mu }\right)\hfill \\ \hfill & +\frac{C_5^V}{m_N^2}\left(g^{\beta \mu }q·k-q^{\beta }{k}^{\mu }+C_6^V{g}^{\beta \mu }\right)]{\gamma }_5,\hfill \\ \hfill {\left(j_a^{\mu }\right)}_A& ={\left(k_{\pi }^{\prime }\right)}^{\alpha }{G}_{\alpha \beta }\left(p_{\mathsf{\Delta }}\right)\left[\frac{C_3^A}{m_N}\left(g^{\beta \mu }\cancel{q}-q^{\beta }{\gamma }^{\mu }\right)+\frac{C_4^A}{m_N^2}\left(g^{\beta \mu }q·p_{\mathsf{\Delta }}-q^{\beta }{p}_{\mathsf{\Delta }}^{\mu }\right)+C_5^A{g}^{\beta \mu }+\frac{C_6^A}{m_N^2}{q}^{\mu }{q}^{\alpha }\right],\hfill \end{array}$$

(44)

where k is the momentum of the initial nucleon which absorbs the incoming momentum $\tilde{q}$ and $p_{\mathsf{\Delta }}=\tilde{q}+k$, yielding $p_{\mathsf{\Delta }}^0=e\left(\mathbf{k}\right)+\tilde{\omega }$. We introduced $\tilde{q}=(\tilde{\omega },\mathbf{q})$ to account for the fact that the initial nucleons are off-shell. A similar definition can be written down for $j_b^{\mu }$; more details are reported in Refs. [35,107]. For $C_3^V$, we adopted the model of Ref. [108], yielding

$$\begin{array}{c}\hfill C_3^V=\frac{2.13}{{(1-q^2/M_V^2)}^2}\frac{1}{1-q^2/\left(4M_V^2\right)},\end{array}$$

(45)

with $M_V=0.84$ GeV. Following the discussion of Ref. [105], 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 Equation (44), which is the term proportional to $C_5^A$ defined as [109]

$$\begin{array}{c}\hfill C_5^A=\frac{1.2}{{(1-q^2/M_{A\mathsf{\Delta }})}^2}\times \frac{1}{{(1-q^2/\left(3M_{A\mathsf{\Delta }}\right))}^2},\end{array}$$

(46)

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

The Rarita–Schwinger propagator

$$\begin{array}{c}\hfill G^{\alpha \beta }\left(p_{\mathsf{\Delta }}\right)=\frac{P^{\alpha \beta }\left(p_{\mathsf{\Delta }}\right)}{p_{\mathsf{\Delta }}^2-M_{\mathsf{\Delta }}^2},\end{array}$$

(47)

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

$$\begin{array}{cc}\hfill \mathrm{\Gamma }\left(p_{\mathsf{\Delta }}\right)& =\frac{{\left(4f_{\pi N\mathsf{\Delta }}\right)}^2}{12\pi m_{\pi }^2}\frac{{\left|\mathbf{d}\right|}^3}{\sqrt{s}}(m_N+E_d)R\left({\mathbf{r}}^2\right).\hfill \end{array}$$

(48)

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

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

(49)

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

$$R\left({\mathbf{r}}^2\right)=\left(\frac{{\mathrm{\Lambda }}_R^2}{{\mathrm{\Lambda }}_R^2-{\mathbf{r}}^2}\right),$$

(50)

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 ${\mathrm{\Lambda }}_R^2=0.95m_N^2$, is introduced to improve the description of the experimental phase shift ${\delta }_{33}$ [85]. The medium effects on the $\mathsf{\Delta }$ propagator are accounted for by modifying the decay width as

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathsf{\Delta }}\left(p_{\mathsf{\Delta }}\right)\to {\mathrm{\Gamma }}_{\mathsf{\Delta }}\left(p_{\mathsf{\Delta }}\right)-2\mathrm{Im}\left[U_{\mathsf{\Delta }}(p_{\mathsf{\Delta }},\rho ={\rho }_0)\right],\end{array}$$

(51)

where $U_{\mathsf{\Delta }}$ is a density-dependent potential obtained from a Bruckner–Hartree–Fock calculation using a coupled-channel $NN\oplus N\mathsf{\Delta }\oplus \pi NN$ model [111,112,113,114], and we fixed the density at the nuclear saturation value ${\rho }_0$ = 0.16 fm${}^3$. For a detailed analysis of medium effects in the MEC contribution for electron–nucleus scattering, see Ref. [107]. One key point to be made from the point of inclusive and even semi-exclusive observables is that the one- and two-body currents contribute coherently, i.e., their interference terms are nonzero. The interference between one- and two-body currents leading to two-nucleon emission has been found to be small [103], but the same interference also contributes to single-nucleon final states [115,116]. The impact of the latter interference on the SF formalism remains to be studied.

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

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

5. Spectral Function

In the factorization scheme of Section 4, the spectral function is the central object containing all the dynamical information about the nucleus. The spectral function of a nucleon with isospin ${\tau }_k=p,n$ and momentum $\mathbf{k}$ can be written as

$$\begin{array}{cc}\hfill P_{{\tau }_k}(\mathbf{k},E)& =\sum _n|\langle {\mathrm{\Psi }}_0|\left[\right|k\rangle \otimes |{\mathrm{\Psi }}_n^{A-1}{\rangle \left]\right|}^2\delta (E+E_0-E_n^{A-1}).\hfill \end{array}$$

(52)

where E is the excitation energy of the remnant, $|k\rangle$ is the single-nucleon state, $|{\mathrm{\Psi }}_0\rangle$ is the ground state of the Hamiltonian in Equation (5) with energy $E_0$, while $|{\mathrm{\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 E

$$n_{{\tau }_k}\left(k\right)=\int dE{P}_{{\tau }_k}(\mathbf{k},E),$$

(53)

and the proton and neutron spectral functions are normalized so that

$$\begin{array}{cc}\hfill & \int dE\frac{d^{3}k}{{\left(2\pi \right)}^3}{P}_p(\mathbf{k},E)=Z,\hfill \\ \hfill & \int dE\frac{d^{3}k}{{\left(2\pi \right)}^3}{P}_n(\mathbf{k},E)=A-Z.\hfill \end{array}$$

(54)

We can rewrite the spectral function as a sum of a mean field (MF) and a correlation (corr) term. The MF piece contains the shell structure with nucleons occupying orbitals obeying the Pauli principle and contributes to the low k and E region. On the other hand, the correlation term comes from pairs and triplets of interacting nucleons with low center-of-mass momentum but large relative momentum above $k_f$. A large body of experimental evidence from $(e,e^{\prime }p)$ data has shown that the correlation piece leads to a depletion of the single-nucleon strength in the MF region by approximately 20% and is essentially nucleus-independent [3,21,117,118,119,120].

Many calculations of the spectral function for finite nuclei are available from a combination of fits to $(e,e^{\prime }p)$ cross sections and theoretical calculations. The spectral function of Benhar et al. obtains the mean field piece from fits to exclusive electron scattering data, and computes the correlation piece from CBF theory for nuclear matter [121,122]. The local density approximation (LDA) is used to extrapolate the correlation piece to finite nuclei by convoluting the correlation component of the nuclear spectral function $P_{\mathrm{NM}}^{\mathrm{corr}}$ with the density profile of the nucleus ${\rho }_A\left(\mathbf{R}\right)$ [123]; it reads

$$P_{\mathrm{CBF}}^{\mathrm{corr}}(\mathbf{k},E)=\int d^{3}R{\rho }_A\left(\mathbf{R}\right)P_{\mathrm{NM}}^{\mathrm{corr}}(\mathbf{k},E;{\rho }_A\left(\mathbf{R}\right)).$$

(55)

In addition to the CBF, the spectral function of nuclear matter and finite nuclei has been computed within the Self-Consistent Green’s Function approach. The latter is a so-called abinitio method that starts from a nuclear Hamiltonian such as Equation (5) with NNLO chiral interactions [124,125]. The SCGF method involves an iterative calculation of the Green’s function’s imaginary component, which is directly related to the one-body spectral function. This technique can be extended to open-shell nuclei and has polynomial scaling with the number of particles, making it feasible for systems with up to $A=100$ [126]. Both CBF and SCGF spectral functions have been used to compute inclusive electron and neutrino scattering cross sections, and have been shown to provide good agreement with electron data when final-state interactions are taken into account [127]. Even though the two spectral functions are obtained from different nuclear interactions, these calculations show that the two many-body approaches produce similar results for inclusive cross sections. Exclusive predictions will most likely be necessary to distinguish the two models.

In this work, we focus on a novel Quantum Monte Carlo (QMC) calculation of the one- and two-body spectral functions for ${}^{12}$C. We begin with the MF piece of the one-body SF for the ${\tau }_k=p$ case. The MF contribution is obtained by considering only the bound $A-1$ states of the remnant nucleus

$$\begin{array}{cc}\hfill P_p^{\mathrm{M}F}(\mathbf{k},E)& =\sum _n|\langle {\mathrm{\Psi }}_0|\left[\right|k\rangle \otimes |{\mathrm{\Psi }}_n^{A-1}{\rangle \left]\right|}^2\hfill \\ \hfill & \times \delta \left(E-B_0^A+B_n^{A-1}-\frac{k^2}{2m_{A-1}}\right),\hfill \end{array}$$

(56)

where $B_0^A$ and $B_n^{A-1}$ are the binding energies of the initial and the bound $A-1$ spectator nucleus with mass $m_{A-1}$. The momentum-space overlaps ${\mathrm{\Psi }}_0|\left[\right|k\rangle \otimes |{\mathrm{\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 following transitions [128,129]:

$$\begin{array}{cc}\hfill {}^{12}\mathrm{C}\left(0^{+}\right)& {\to }^{11}\mathrm{B}(3/2^{-})+p\hfill \\ \hfill {}^{12}\mathrm{C}\left(0^{+}\right)& {\to }^{11}\mathrm{B}(1/2^{-})+p\hfill \\ \hfill {}^{12}\mathrm{C}\left(0^{+}\right)& {\to }^{11}\mathrm{B}{(3/2^{-})}^{*}+p.\hfill \end{array}$$

The calculation of the s-shell mean-field contribution involves nontrivial difficulties for the VMC method, as it would require an evaluation of the spectroscopic overlaps for the transitions to all the possible excited states of ${}^{11}\mathrm{B}$ with $J^P=(1/2^{+})$. To overcome this limitation, we used the VMC overlap associated with the ${}^4\mathrm{He}\left(0^{+}\right)\to {}^3\mathrm{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 Refs. [128,129].

The correlation contribution to the SF is given by

$$\begin{array}{cc}\hfill P_p^{\mathrm{c}orr}& (\mathbf{k},E)=\sum _n\int \frac{d^3{k}^{\prime }}{{\left(2\pi \right)}^3}|\langle {\mathrm{\Psi }}_0\left|\left[\right|k\rangle \right|k^{\prime }\rangle |{\mathrm{\Psi }}_n^{A-2}{\rangle \left]\right|}^2\hfill \\ \hfill & \times \delta (E+E_0-e\left({\mathbf{k}}^{\prime }\right)-E_n^{A-2})\hfill \\ \hfill & ={\mathcal{N}}_p\sum _{{\tau }_{k^{\prime }}=p,n}\int \frac{d^3{k}^{\prime }}{{\left(2\pi \right)}^3}[n_{p,{\tau }_{k^{\prime }}}(\mathbf{k},{\mathbf{k}}^{\prime })\hfill \\ \hfill & \times \delta \left(E-B_0-e\left({\mathbf{k}}^{\prime }\right)+{\overline{B}}_{A-2}-\frac{{(\mathbf{k}+{\mathbf{k}}^{\prime })}^2}{2m_{A-2}}\right)],\hfill \end{array}$$

(57)

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 ${\overline{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. [130], 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 were correctly recovered. Below in Figure 5, we show the ${}^{12}$C single-nucleon momentum distribution derived using the above procedure. Figure 5 shows the effect of the different prescriptions for calculating the s-wave overlaps, with the above prescription resulting in an increased normalization of the SF compared to the harmonic oscillator and Wood–Saxon potentials. In Figure 6, we directly compare the QMC and CBF ${}^{12}$C spectral functions by comparing their one-dimensional momentum and removal energy distributions. While the two SFs have very similar removal energy distributions, their momentum distributions show distinct behavior at small and large nucleon momenta. Although these discrepancies only cause minor variations in the inclusive cross section, it is anticipated that they will be more significant in exclusive cross sections where the outgoing nucleon is measured. This will be explored in future studies.

For the contribution of multi-nucleon currents to the cross section, a two-nucleon spectral function is needed. As mentioned previously, in infinite nuclear matter a factorization of the two-nucleon momentum distribution into the product of two single-nucleon momentum distributions can be made. This factorization assumes no long-range correlations are present and throws away correlations between the two struck particles. We go beyond this approximation by explicitly by using the two-nucleon spectral momentum distribution to build the two-nucleon spectral function. We include only the mean field contribution, e.g., we neglect contributions where more than two nucleons are emitted, which reads

$$\begin{array}{cc}\hfill P_{{\tau }_k,{\tau }_k^{\prime }}^{\mathrm{M}F}(\mathbf{k},{\mathbf{k}}^{\prime },E)& =n_{{\tau }_k,{\tau }_{k^{\prime }}}(\mathbf{k},{\mathbf{k}}^{\prime })\hfill \\ \hfill & \times \delta \left(E-B_0+{\overline{B}}_{A-2}-\frac{{\mathbf{K}}^2}{2m_{A-2}}\right),\hfill \end{array}$$

(58)

where $\mathbf{K}=\mathbf{k}+{\mathbf{k}}^{\prime }$ is the total momentum of the pair.

6. GFMC and SF Comparisons

Recently, the authors of Ref. [37] computed flux-folded differential cross sections for MiniBooNE and T2K experiments [131,132] using the QMC spectral function outlined above. Under control systematics, the calculation of the two-body contribution is required for disentangling the effect of the axial form factor. The QMC spectral function can be directly compared with GFMC predictions because it is derived from the same underlying Hamiltonian and currents. Comparisons with experiments have shown that the predictions are consistent with the data and show tension between the results obtained from adopting the LQCD and phenomenological form factors displayed in Figure 1 [75]. The results for selected angular bins for T2K kinematics are shown in Figure 7. The shown uncertainty bands propagate the uncertainty on the axial nucleon form factor derived using the z expansion, where the different coefficients have been fitted to the deuterium bubble chamber data of Ref. [69]. In the case of the GFMC, the uncertainty coming from the inversion of Euclidean responses is also included. Both approaches provide a similar description of the data, albeit the contribution of two-body currents peaks is shifted in the two approaches. This can be ascribed to different motivations. Firstly, the SF results include explicitly the contribution of $\mathsf{\Delta }$ excitations in the two-nucleon knockout process, leading to a peak at smaller lepton momenta, while the GFMC results use a static $\mathsf{\Delta }$ treatment. Secondly, the GFMC results also account for the interference between the two-body and one-body currents, which would also lead to an enhancement in the vicinity of the quasi-elastic peak. Such an enhancement is clearly seen in Figure 7, and in the electromagnetic and electroweak responses [27,28]. While these observations support the one- and two-body current interference, it is impossible to disentangle this contribution directly in the GFMC results. The calculations in the nuclear matter and relativistic mean field calculations of Refs. [115,116] also find that the transverse enhancement observed in electron scattering is primarily due to the constructive interference between one- and two-body currents, leading to single-nucleon knockout final states.

Recently [34], relativistic corrections to GFMC calculations for flux-averaged neutrino cross sections have been determined using the method described in Section 3.1.1. The influence on the T2K results shown in Figure 7 is small and generally falls within the uncertainty bands due to the axial form factor. For MINER$\nu$A data [38] taken with the medium-energy NuMI beam, which peaks at around $6\mathrm{GeV}$ [133], relativistic corrections are crucial. The charged-current flux-averaged cross section is presented in terms of muon momentum parallel and perpendicular to the beam direction

$$p_{\Vert }=\left|{\mathbf{p}}_{\mu }\right|cos{\theta }_{\mu },$$

(59)

and

$$p_{\perp }=\left|{\mathbf{p}}_{\mu }\right|sin{\theta }_{\mu }=\sqrt{{\mathbf{p}}_{\mu }^2-p_{\Vert }^2},$$

(60)

respectively.

In this comparison, the routinely used dipole parameterization of the axial form factor with ${\mathrm{\Lambda }}_A=1.03$ GeV has been used. For the GFMC results, the error band includes statistical errors combined with the error from the inversion of Euclidean responses. The SF results do not include an error estimate. Relativistic corrections are included in the GFMC results by performing the calculation in the active-nucleon Breit frame (ANB) as discussed in Section 3.1.1. The incorporation of relativistic effects leads to a nearly halved cross section for low $p_{\Vert }$, with the discrepancy gradually decreasing as $p_{\Vert }$ increases. We observe that the momentum transfer is constrained such that $q>p_{\perp }$, and smaller $p_{\Vert }$ bins generally permit higher energy and consequently larger q contributions at small $p_{\perp }$, thus explaining this behavior. The emergence of high-$p_{\perp }$ (i.e., high-q) tails can be understood as the response’s narrowing in terms of energy transfer compared to nonrelativistic outcomes, resulting in strength redistribution within the available phase space at large q.

Given the inclusion of large q values in the MINER$\nu$A calculations and the substantial impact of relativistic corrections, a consistency check is warranted. In Ref. [34], the GFMC results for MINER$\nu$A kinematics obtained including only the one-body current contribution have been compared to other approaches that are either manifestly relativistic [134] or include relativistic corrections [13,40,135,136] and found to agree with the theoretical curves. Here, Figure 8 compares the GFMC results to the SF calculations including both the one- and two-body contributions in Figure 8. The agreement between the one-body contribution in the GFMC and SF approaches is evident when the former are computed in the ANB. The total increase in the cross section due to two-body contributions is twice as large in the SF calculations compared to the GFMC. This difference can be attributed to the same motivations discussed for the T2K results.

Lastly, we note that the GFMC nonrelativistic calculations exhibit better conformity with experimental data compared to those incorporating relativistic effects. However, considering the energy distribution of the medium-energy NuMI beam in the MINER$\nu$A experiment, it is expected that contributions beyond quasi-elastic scattering are significant, even when events with detectable mesons are excluded from the experimental analysis. Specifically, there are instances where pions produced at the interaction vertex are either absorbed or remain undetected. Thus, theoretical calculations that neglect pion-production mechanisms should yield results lower than experimental data. This aligns with the case where relativistic effects are considered, while their omission leads to un-physically large cross sections.

In this review, we also consider inclusive electron scattering data on ${}^{12}$C in Figure 9, which allows one to disentangle the different energy regions more clearly. The two kinematics under consideration have been deliberately selected to include only responses with $q<700$ MeV. These specific values align with the range for which the GFMC responses have already been computed. Similarly to the neutrino case, in the GFMC calculations two-body currents provide an enhancement in the quasi-elastic region. A comparison between the LAB frame results using purely relativistic kinematics (depicted by the blue dotted curve) and the ANB curve (solid blue) reveals important insights. Relativistic corrections cause a shift of the peak towards smaller $\omega$ values, a reduction in width and an increase in the height of the quasi-elastic peak. The one-body contribution computed in the ANB frame, displayed by a solid red line, agrees fairly well with the SF one-body contribution displayed in the upper panels. Overall agreement with the data improves by including relativistic corrections to the GFMC results. However, it is worth noting that the absence of $\pi$-production contributions makes it difficult to draw definitive conclusions without considering that term.

The static treatment of the $\mathsf{\Delta }$ propagator restricts the significance of two-body currents in the “dip” region, located between the quasi-elastic and pion-production peaks. Incorporating explicit dynamical degrees of freedom in GFMC calculations is more challenging, particularly in terms of evaluating the Euclidean responses while fixing the current operator’s dependence at multiple values of $\omega$.

The total QMC SF results encompass the incoherent sum of one-nucleon and two-nucleon contributions. We include in the calculation the effect of final-state interactions by convoluting the computed cross sections with a folding function which both shifts and redistributes strength from the peak to the tails [136]. Two-body currents give a minor enhancement in the quasi-elastic peak region, but a strong enhancement in the “dip” region. Additionally, we include CBF SF results for the one-body cross section for comparison, which show similar trends to the QMC one-body cross section as expected. However, the QMC SF result notably under-predicts the data in the region of the quasi-elastic peak at $E_{\mathrm{beam}}=620 \mathrm{MeV}$. Investigating the interference between one -and two-body currents and its impact on these results will be a subject of future investigation.

As a general remark, one can choose to apply either the GFMC or the spectral function approach depending on the kinematics and process under investigation. The selection depends on the specific requirements of the study. However, it is important to ensure that the results obtained from both methods are consistent in the transition regions where both approaches are expected to work.

7. Conclusions

Neutrino oscillation experiments cover a broad range of energies, from a few MeV to tens of GeV, where different reaction mechanisms involving various degrees of freedom (nucleons, pions, quarks, etc.) are active. Microscopic approaches such as Green’s Function Monte Carlo (GFMC) and Coupled Cluster have been successful in describing lepton–nucleus cross sections in the MeV energy region [140,141,142]. However, to address the higher energies relevant for DUNE and include explicit pion degrees of freedom, different methods relying on a factorization of the hadronic final state, such as the Spectral Function (SF), the Short-Time Approximation [143,144] and the Relativistic Mean Field approach [12,134], have proven successful in reproducing electron scattering data for different kinematics.

Providing a realistic estimate of the theoretical uncertainty of the prediction in the neutrino–nucleus cross section, which must be propagated in the extraction of neutrino oscillation parameters, requires assessing the error associated with the input used in the calculations and with the many-body method used. In this review, we highlight that different choices can be made to define the nuclear forces adopted to describe the wave function of the target and remnant nuclei, either using semi-phenomenological approaches or chiral effective field theories. Following the choice of the nuclear forces, different current operators can also be constructed. Another source of uncertainty is connected to the form factors entering these currents. In Ref. [37] a study of the dependence of the neutrino–nucleus cross section results from the axial form factor adopted in the one-body current operator has been carried out using the GFMC and the SF approaches and a tension between the results obtained for the LQCD and phenomenological form factors has been observed. The results of Ref. [37] indicate that while significant progress has been made in the determination of the axial and vector form factors entering in the one-body current operator, more work will be required in the future for the determination of the form factors entering the two-body currents, particularly for those contributions with $\mathsf{\Delta }$ degrees of freedom.

A two-fold strategy can be employed to comprehend the error associated with using a factorization scheme in the spectral function approach and nonrelativistic kinematics in GFMC. Firstly, relativistic corrections can be incorporated by working in a reference frame that minimizes them in the GFMC responses [34,140]. Secondly, Quantum Monte Carlo (QMC) techniques can be used to derive one- and two-nucleon spectral functions. Comparing the results obtained from these two approaches can help estimate the error associated with the many-body method. Numerous studies have investigated this comparison. In this review, we present unpublished results that demonstrate electron–carbon cross section comparisons and neutrino–nucleus cross sections for the MINER$\nu$A experiment. In the comparison with MINER$\nu$A Medium-Energy CCQE-like data, the effect of relativistic corrections to the GFMC results are substantial, yielding a quenching of the results up to 50% of the initial strength. We observe a reasonable agreement between the GFMC and QMC SF results. For the electron scattering cross section, we also analyzed the dependence of the results from the many-body method adopted to derive the spectral function; in particular, we compared the QMC and Correlated Basis Function results and found a very good agreement between them. Looking at a fixed energy beam allows one to better separate the contribution of the different reaction mechanisms. In this case, the difference between the two-body contributions obtained within the two approaches is apparent and it has to be attributed to the different treatment of the $\mathsf{\Delta }$ propagator in the GFMC and the lack of one- and two-body current interference in the SF approach. The inclusion of relativistic corrections in the GFMC results leads to better agreement with data. As there is a large amount of electron scattering data in the region of $300<q<700$ MeV, future studies that directly compare the GFMC results with differential electron scattering data for carbon can be performed. A robust method for estimation of the uncertainty able to account for all the different aspect of the calculation is required to match the unprecedented accuracy of neutrino experiments, some preliminary steps toward this direction have been discussed in this review using the GFMC and SF methods.