Ab initio Calculations of Lepton-Nucleus Scattering

The advent of high precision measurements of neutrinos and their oscillations calls for accurate predictions of their interactions with nuclear targets utilized in the detectors. Over the past decade, ab initio methods based on realistic nuclear interactions and current operators were able to provide accurate description of lepton-nucleus scattering processes. Achieving a comprehensive description of the different reaction mechanisms active in the broad range of energies relevant for oscillation experiments required the introduction of controlled approximations of the nuclear many-body models. In this review, we give an overview of recent developments in the description of electroweak interactions within different approaches and discuss the future perspectives to support the experimental effort in this new precision era.


INTRODUCTION
Understanding neutrino properties and interactions is the main focus of the world-wide accelerator-based neutrino-oscillation program. The new generation of short [1] and longbaseline [2] neutrino experiments-such as the Deep Underground Neutrino Experiment (DUNE)-will address fundamental questions and play a key role in the search of physics beyond the Standard Model. In particular, the absolute scale of neutrino masses and the presence of charge parity (CP) violation in the leptonic sector, which may contribute to our understanding of the matter-antimatter asymmetry of the universe, will be determined. The existence of a fourth (sterile) neutrino will be investigated, this could explain the excess of electron neutrinos from charged current quasi-elastic events reported by the MiniBooNe collaboration [3].
Studying neutrino-nucleus interactions in the energy region of interest for oscillation experiments is a multi-scale problem. In fact, the experiments [4][5][6][7][8][9] are sensitive to a broad range of energy where different reaction mechanisms contribute [10,11]. Quasielastic scattering is dominant for energies of hundreds of MeVs, in this case the lepton interacts with individual bound nucleons. In addition to this, there are corrections accounting for coupling of the probe to interacting nucleons, belonging to a correlated pairs on connected via two-body currents.
For larger values of the energy, a struck nucleon can be excited to a baryon resonance state and subsequently decays into pions, or give rise to deep-inelastic scattering (DIS) processes. Constructing a framework able to describe these diverse reaction mechanisms on the same footing is a formidable nuclear-theory challenge. Nuclear effective field theory (EFT) provides a systematic way to construct nuclear interactions and currents performing a low-momentum expansion. In addition, an estimate of the theoretical uncertainty-which will be crucial for the neutrino data analyses-can be properly assessed. In the last decade there has been a tremendous progress in the field of nuclear ab initio methods made possible by the increasing availability of computing resources and the development of new algorithms. Within these approaches the nucleus is treated as a collection of nucleons interacting with each other via two-and three-body forces obtained within nuclear EFTS. The interaction with external electroweak probes is described by one-and two-body effective currents that are consistent with the nuclear potentials [12][13][14][15][16][17].
Among these ab initio many-body methods, the Green's Function Monte Carlo (GFMC) method utilizes quantum Monte Carlo (QMC) techniques [18] to perform first-principle calculations of nuclei up to 12 C [19,20]. More recently the GFMC method has also been applied to the calculation of the electroweak sum rules and response functions of 4 He and 12 C, including one and two-body currents [21][22][23]. In order to do that, integral transform techniques have been utilized. The accuracy of the inversion procedure adopted within the GFMC to obtain the response function from its integral transform has been compared and benchmarked with the Lorentz Integral Transform method [24]. The GFMC electroweak responses computed in the quasielastic sector are virtually exact for low and moderate values of the momentum transfer. Initial and final state correlations are fully retained within this approach [25,26]. In order to extend the predictive power of this approach, relativistic effects in nuclear kinematics were included [24] leading to an excellent agreement with electron scattering data off 4 He. However, within this approach it is not possible to access more exclusive channels, the calculations are fully inclusive. The inclusion of explicit pion degrees of freedom in the nuclear wave function [27], necessary to describe the resonance production region, is extremely complicated and its achievement is still a long way ahead.
The short-time approximation (STA) method has been recently proposed to overcome some of the limitations of GFMC [28]. In particular, this approach allows to compute both the inclusive and exclusive response of nuclei in the highenergy (short-time) limit-corresponding to the Fermi energy and above-utilizing realistic nuclear interactions and currents. In the STA the full ground-state dynamics is retained while the hadronic final state is factorized at the two-nucleon level, this approximation is expected to be valid at high-energy and momentum transfer. It allows to account for the final state interactions of the pairs involved in the interaction vertex and to incorporate two-nucleon correlations and currents as well as their interference, which are known be sizable in the GFMC results [22][23][24][25][26]. While keeping consistently two-body physics and ensuing quantum interference contributions, the STA is not expected to accurately model the correct physics for low-lying excitations or collective behavior like giant resonances. A good agreement with the GFMC electromagnetic responses of 4 He is observed after enforcing the correct threshold behavior. Since the STA involves only two active nucleons, it suitable to be improved and include relativistic currents and kinematics as well pion production channels.
In the past years, the framework based on the impulse approximation (IA) and realistic spectral-functions (SFs) has been largely utilized to describe electron-nucleus scattering data in the limit of moderate and high momentum transfer [29,30]. This scheme combines a realistic description of the initial target state with a fully-relativistic interaction vertex and kinematics. This is achieved factorizing the hadronic final state in terms of a free nucleon state and enclosing all the information on nuclear structure in the SF. The latter does not depend on the momentum transfer and can be computed within nonrelativistic nuclear many-body approaches. In this review we consider two nuclear SFs, derived within the correlated basis function (CBF) formalism [31] and the self-consistent Green's function (SCGF) theory [32,33]. These two approaches utilize different nuclear forces and involve different approximations in each of the SF calculations. Within the factorization scheme lepton-nucleus scattering is rewritten as an incoherent sum of elementary processes involving individual nucleons. This framework has been extended and generalized to include twonucleon emission processes induced by relativistic mesonexchange currents [34] and applied to calculate the electroweak inclusive cross sections of carbon and oxygen [35,36]. In order to tackle the resonance production region, the electroweak pion production amplitudes generated within the dynamical coupled-channel (DCC) model [37][38][39] have been included in the IA formalism. The results obtained convoluting the DCC elementary amplitudes with the CBF SF will also be reported.
In this review, we present lepton-nucleus interaction results obtained within different many-body methods. In section 2 we outline the formalism and define the electroweak cross sections, nuclear interactions and currents. In section 3 we discuss the integral transform techniques utilized in the GFMC to obtain the nuclear response functions. Sections 4 and 5 are devoted to the STA and extended factorization scheme, respectively. We present recent results obtained within each of the aforementioned approaches for different nuclei. In section 6 we state out conclusions and discuss future directions.

FORMALISM
We consider the scattering of an initial electron of fourmomentum k = (E, k) off a nucleus A at rest; in the final state only the outgoing electron with momentum k = (E ′ , k ′ ) is detected. The inclusive double differential cross section for this process can be written in Born Approximation as where α ≃ 1/137 is the fine structure constant and ′ is the scattering solid angle in the direction specified by k ′ . The energy and the momentum transfer are denoted by ω and q, respectively, with Q 2 = −q 2 = q 2 − ω 2 . The lepton tensor is fully determined by the lepton kinematical variables. It can be written neglecting the electron mass as The hadronic tensor describes the transition between the initial and final nuclear states | 0 and | f , with energies E 0 and E f .
For nuclei of spin-zero, it can be written as where we sum over all hadronic final states and J µ (q, ω) is the nuclear current operator whose structure will be discussed in detail in section 2.1. For inclusive processes, the cross section of Equation (1) only depends on the two response functions, R L (q, ω) and R T (q, ω), which describe the interactions with a virtual photon longitudinally (L) and transversely (T) polarized. It is given by where the lepton kinematical factors are given by and is the Mott cross section depending on the scattering angle θ . The extension to neutral-and charge-current electroweak processes is straightforward. Let us consider the scattering of a neutrino (ν ℓ ) or an anti-neutrino (ν ℓ ) off a nuclear target. In analogy with the electromagnetic case, we denote by k = (E, k) and k = (E ′ , k ′ ) the momentum of the initial and outgoing lepton; the double-differential cross section can be written as [40,41] dσ where G = G F and G = G F cos θ c for neutral current (NC) and charge current (CC) processes, respectively, with cos θ c = 0.97425 [42] and for the Fermi coupling constant we use G F = 1.1803 × 10 −5 [43]. The leptonic tensor contains an extracontribution proportional to the Levi Civita where the + (−) sign is for ν ℓ (ν ℓ ) initiated reactions. The hadronic tensor is the same as Equation (3) but the current operator will have a vector and axial component, its explicit expression will be given in section 2.1. We perform the Lorentz contraction of the leptonic and hadronic tensor of Equation (7), yielding Taking the three-momentum transfer along the z axis and the total three-momentum in the x − z plane, we can define q and Q as and write the leptonic kinematical factors aŝ with m 2 ℓ = k ′ 2 the mass of the outgoing lepton. The five electroweak response functions correspond to different components of the hadro tensor Note that electron and neutrino scattering cross sections are written in a similar fashion as a contraction of the leptonic and the hadronic tensor. This analogy will become even more apparent in section 2.1 where we introduce the explicit expression of the current operators and use the conserved vector current (CVC) hypothesis to connect the vector part of the electroweak current with the electromagnetic one. For this reason, any model of neutrino-nucleus scattering has to be capable of describing electron-scattering cross sections first [44].

Nuclear Hamiltonian and Current Operator
The internal structure of nuclei and their reactions can be described within non-relativistic many body approaches utilizing an Hamiltonian in which the nucleons are the only active degrees of freedom. Its general expression is given by where p i is the momentum of the i-th nucleon having mass m N , while the potentials v ij and V ijk model the nucleon-nucleon (NN) and three-nucleon (3N) interactions, respectively. A model of NN interaction has to be constrained by the large number of NN scattering data available. Currently, very accurate results have been obtained from both phenomenological approaches and chiral effective field theory able to accurately fit these data. The Argonne v18 is a finite, coordinate space NN potential that has been fit to the full Nijmegen phase-shift database and to low energy scattering parameters and deuteron properties.
It can be written as a product of radial functions and spinisospin operator v ij = p=1,...,18 v p (r ij )O p ij (14) where r ij = |r i −r j |. The first 14 operators are charge independent (corresponding to the older Argonne v14 model) where σ and τ are the Pauli matrices operating over the nucleon spin and isospin degrees of freedom, respectively, is the relative angular momentum, and S ij = 1 2 (σ i + σ j ) is the total spin. The terms p = 15 . . . 18 are included in the formulation of the Argonne potential to account for small violations of charge symmetry [45]. In addition to the NN interactions, also phenomenological 3N interactions have been developed. In particular the Illinois 7 (IL7) interaction is expressed as a sum of a two-pion-exchange P-wave term (Fujita-Miyazawa), a two-pion-exchange S-wave contribution, a three-pion-exchange contribution, and a 3N central interaction [46]. These phenomenological NN and 3N interactions were successfully utilized in QMC calculations of nuclear properties, neutron drops, and neutron-star matter [18,[47][48][49][50].
In the last decades there has been a tremendous progress in development of chiral EFT interactions as proven by the availability of a number of different potentials [13,[51][52][53][54]. However, it is only recently with their formulation in coordinate space that the application of chiral forces within QMC approaches became possible. These new local potentials have been computed up to N 2 LO with consistent 3N terms in Gezerlis et al. [55,56], the -isobar degrees of freedom were fully accounted for in the derivation of Piarulli et al. [57,58].
In analogy with the nuclear Hamiltonian, the current operator can be expressed as a sum of one-and two-body terms The one-body electromagnetic current is given by where p and p ′ are the initial and final nucleon momentum. The isoscalar (S) and isovector (V) form factors, F 1 and F 2 , are given by combination of the Dirac and Pauli ones, F 1 and F 2 , as where τ z is the is the isospin operator and The Dirac and Pauli form factors can be expressed in terms of the electric and magnetic form factors of the proton and neutron as with τ = Q 2 /4m 2 N . Therefore, the electromagnetic current can be schematically written as J µ EM = J µ γ ,S + J µ γ ,z where the first is the isoscar term and the second is the isovector multiplied by the isospin operators τ z . The one-body charge and current operator are obtained from the non-relativistic reduction of the covariant operator of Equation (17) including all the terms up to 1/m 2 N in the expansion. It leads to the following expression for isoscalar charge, transverse (⊥) and longitudinal ( ) to q component of the current operator Note that the last relation has been obtained from current conservation relation discussed in Equation (30). Analogously to Equation (19), the isoscalar and isovector component of the electric and magnetic form factors are written as The isovector contribution to the current J µ γ ,z is obtained by replacing G S E,M → G V E,M τ z . The electroweak interactions of a neutrino or anti-neutrino with the hadronic target can be induced by both CC and NC transitions. The current operator is written as a different combination of vector and axial terms. Note that in both cases the Conserved Vector Current (CVC) hypothesis allows to relate the vector form factor to the electromagnetic ones where θ W is the Weinberg angle (sin 2 θ W = 0.2312 [42]). The fully relativistic expression of the axial one-body current operator reads For CC processes the axial and pseudo-scalar form factors can be written as with τ i,± = (τ i,x ± τ i,y )/2 being the isospin raising-lowering operator.
The dipole parametrization for the axial form factor is routinely adopted in most of the calculations, it reads where the nucleon axial-vector coupling constant is assumed to be g A = 1.2694. The impact of the Q 2 dependence of the axial form factor on neutrino-nucleus cross-section predictions has been discussed in Aguilar-Arevalo et al. [59] and Bernard et al. [60]. The validity of the dipole parametrization has been questioned and in these regards both lattice-QCD calculations [61] and the so called "z-expansion" analysis [62] has been recently proposed. The pseudo-scalar form factor is obtained by using Partially Conserved Axial Current (PCAC) arguments to write it in terms of the axial one For NC transitions, we report the non-relativistic reduction of the charge-and axial-current operator [63] (for brevity we neglect order 1/m 2 N terms) The CC non-relativistic current is obtained by substituting in the non-relativistic isovector term j µ γ ,z and j µ a,z of Equation (28), τ i,z /2 → τ i,± and adding a pseudoscalar contribution The electromagnetic current must satisfy the continuity equation which links the two-body exchange current operator to the NN interaction and leads to separate continuity equations for the one-and two-body current operators. Chiral EFT formulations allow to construct electroweak currents which are fully consistent with the nuclear forces entering the hamiltonian. The gauge invariance of the theory implies that the nuclear current operators satisfy the continuity equation with the potentials order by order of the chiral expansion. The calculation of the twobody electromagnetic currents has been the subject of extensive study carried out by different groups [64][65][66][67][68]. Another important advantage of chiral EFT consists on having connected threenucleon interaction and the two-nucleon axial current; their derivation up to one-loop is reported in Krebs et al. [69] and Baroni et al. [70] while the -full expression obtained consistently with the nuclear forces of Piarulli et al. [57,58] can be found in Baroni et al. [71]. The majority of the results that will be presented in this review have been obtained utilizing the aforementioned AV18 potential and related semi-phenomenological currents. The spin and isospin dependence of the NN potential leads to a nonvanishing commutator with the non-relativistic one-body charge operator. In order to satisfy the continuity equation, the socalled "model-independent" two-body currents j MI (ij) have to be introduced; their longitudinal part is directly constrained by the NN interaction while their transverse components are not uniquely defined. The main contributions comes from the onepion and one-rho exchange current operator, their expression is well-known and reported in Rocco et al. [72], Dekker et al. [36] and Schiavilla et al. [73] both in their relativistic and non-relativistic formulation.
In addition to the model-independent two body currents a significant contribution comes from model dependent terms in which the exchange of a pion is followed by the excitation of a -resonance in the intermediate state. Due to the transverse nature of this current, its expression is model dependent. In fact, the form of the vector part is not determined from currentconservation constraints [74].
Analogously to the one-body case, the two-body CC operator is the sum of a vector and axial component. The relativistic expression of the two-body current operator we present in this work has been obtained following the parametrization of Ruiz Simo et al. [75] where the pion-production amplitudes of Hernandez et al. [76] are coupled to a second nucleonic line. Starting from the vector component of j µ CC , the electromagnetic current operators can be obtained using the CVC hypothesis. Adopting the momentum variables specified in Figure 1, the expression of the CC -current reads where k 2 = p ′ − k ′ is the momentum of the π exchanged in the two depicted diagrams, f * = 2.14 and with π N = 1150 MeV and π = 1300 MeV. In the above equations, (k) describes the pion propagation and absorption, while the πNN and πN couplings depend on the form factors F π N (k) and F π NN (k), respectively, accounting for the offshellness of the π and . The operator (I V ) ± = (τ (1) × τ (2) ) ± with ± → x ± iy raises-lowers the isospin components. In Equation (31), j µ a and j µ b denote the N → transition vertices of the left and right diagrams, respectively. They are expressed as where k is the momentum of the initial nucleon, q is the momentum transfer and p = q + k, yielding p 0 = e(k) + ω.
For the left diagram we have where p is the outgoing nucleon four-momentum and p = p − q. The vector and axial form factors, denoted by C V 3 and C A 5 , are obtained from general principles and experimental results as discussed in Hernandez et al. [76] where their explicit expressions is also reported. In the above equations all nucleons are on the mass-shell with the time component p 0 = m 2 N + p 2 . For the -propagator we adopted the Rarita-Schwinger convention where G αβ (p ) = P αβ (p )/(p 2 − M 2 ) is proportional to the spin 3/2 projection operator In order to account for the possible decay of the into a physical πN we replace its real mass M = 1,232 MeV entering the denominator of the free propagator, i.e., p 2 − M 2 , by M − iŴ(p )/2 [72,77]. The decay width Ŵ(p )/2 is not a constant but explicitly depends upon the energy, its expression can be found in Dekker et al. [72].

INTEGRAL TRANSFORM TECHNIQUES
Evaluating the hadronic tensor of Equation (3) is highly nontrivial as it requires to sum over the entire excitation spectrum of the nucleus and to include one-and two-body current operators. Integral transform techniques are extremely useful in reducing the problem to a ground-state one instead of explicitly evaluating each transition amplitude | 0 → | f . We consider the convolution of the response function with a smooth kernel using the closure property f | f f | = 1 a generalized sum rule depending on a continuous parameter σ can be obtained With an appropriate choice of the kernel K, the right-hand side of the above equation can be accurately computed within different ab-initio methods. In order to retrieve the energy dependence of the response function, the integral transform has to be accurately inverted.

Lorentz Integral Transform Technique
The kernel used to compute the integral transform is a Lorentzian where σ is a complex parameter, generally defined for convenience as For simplicity we leave the Lorentz indices implicit and rewrite Equation (39) substituting the Lorentz kernel as implying that the Lorentz Integral Transform (LIT) of the response function is given by the norm of the |˜ state. This state is determined by solving a Schrödinger-like equation for a quantum mechanical bound system for different values of σ I and σ R . In order to obtain the response function, one first computesĒ(q, σ I , σ R ) and then inverts it numerically. This second step is highly non-trivial and belongs to the class of so-called "ill-posed" problems. The devised method consists on: (i) making an ansatz to write down the expression of the response function as a function of a set of parameters c i (ii) the corresponding LIT is computed using this parametrization of the nuclear response function (iii) the parameters c i are determined from a least-squares fit of the LIT computed in (ii) with the one calculated at the beginning,Ē(q, σ I , σ R ). A more detailed discussion of the inversion procedures utilized can be found in Efros et al. [78] and Reiss et al. [79]. A very accurate determination the electromagnetic responses of light nuclei has been achieved, for different values of q, combining the LIT method with the hyperspherical harmonic (HH) formalism [80][81][82]. This expansion method has been successfully utilized to study nuclear, atomic and molecular fewbody systems; the wave function of the system is expanded in a series of products of HH basis functions and hyperradial basis functions allowing for a correct description of the large distance components. The predictive power of the HH approach is limited to relatively light mass number, to overcome this limitation the Couple Cluster (CC) theory has been recently combined with the LIT approach to tackle medium and large mass nuclei [83,84]. The photoabsorption cross sections of 16,22 O and 40 Ca and electromagnetic sum rules have been recently computed using the LIT-CC approach in Simonis et al. [85].

Green's Function Monte Carlo
The Green's Function Monte Carlo (GFMC) is an ab-initio method that allows to predict with high accuracy the structure and low-energy transitions of A ≤ 12 nuclei [19]. This method is utilized to project out the ground state starting from a trial wave where E 0 is a parameter used to control the normalization and τ is the imaginary time. Recently, exploiting integral transform techniques accurate predictions for the electromagnetic response functions of 4 He and 12 C in the quasielastic sector have been obtained. The inclusion of two-body currents yields to a very good agreement between GFMC predictions and experimental data [23,25]. The integral transform of Equation (38) is evaluated using a Laplace kernel and denoted as Euclidean response. Its inelastic contribution is obtained as where ω el is the energy of the recoiling ground state. Using the closure property, the sum over the final states can be removed, the inelastic Euclidean responses can be written as the following ground-state expectation value where F α (q) is the longitudinal elastic form factor and the nucleon form factors entering in the current operator are evaluated at the quasielastic peak Q 2 qe = q 2 − ω 2 qe . In order to obtain the response function the Laplace transform has to be inverted. This is achieved exploiting maximum entropy techniques, as described in Lovato et al. [25] and in the supplemental material of Lovato et al. [86]. The LIT-HH and GFMC results for the longitudinal electromagnetic responses of 4 He have been recently compared to benchmark the inversion procedure used in the two approaches. The calculations presented are based on the AV18 and IL7 combination of twoand three-nucleon potential for the GFMC method and AV18 and UIX for the LIT-HH [45,87]. The results are displayed in Figure 2 for |q| = 300 and 500 MeV, the (red) dashed and solid (blue) lines corresponding to the LIT-HH and GFMC predictions, respectively, are found in good agreement and correctly reproduce the experimental data taken from Carlson et al. [88]. The small discrepancies between the two curves have been discussed in Rocco et al. [24] and they do not originate from the inversion techniques utilized in the two approaches.
Because of the non-relativistic nature of the GFMC approach, it can be safely applied to compute electroweak responses in the high momentum transfer region relevant for neutrino-nucleus scattering. Relativistic corrections are in both one-and two-body current operators up to order 1/m 2 N as discussed in Carlson and Schiavilla [63] for the one-body case. However, the quantum mechanical framework is non-relativistic; an attempt of including relativistic corrections in the kinematics has been recently carried out in Rocco et al. [24]. The strategy consists on performing the calculation in a reference frame that minimizes nucleon momenta in the final state.
In electron-nucleus scattering processes, the quasielastic region is dominated by a one-nucleon knock-out and this condition is satisfied by the active nucleon Breit (ANB) frame in which the target nucleus has a momentum −A q/2. In this reference frame the momentum of the nucleons in the initial state is about −q/2 while the one of the emitted particle is ≃ q/2. The momentum of the final state is higher in all the other reference frames like for example the laboratory (LAB) system where the emitted knocked-out nucleon has a momentum of about q. Therefore, the ANB system can be utilized to minimize relativistic corrections.
In order to compare with experimental data measured in the LAB system, the results from the ANB need to be boosted back to the LAB frame. The solid blue and red curves in Figure 3 performing the calculation of the longitudinal electromagnetic response of 4 He at |q| = 700 MeV in the LAB and ANB frame, respectively, and boosting back to the LAB frame using Lorentz transformation. The two curves peak in different positions and have different strengths, this frame dependence of the results is a clear indication that relativistic effects are sizable for these values of q.
Relativistic effects in the kinematics can be included employing the two-fragment model of Efros et al. [89] which presents strong analogies with the procedure used to determine NN potentials. In this case, the relative momentum of the twonucleon system p 12 is determined in a relativistic fashion and utilized to solve the Schrödinger equation with a non-relativistic kinetic energy E 12 = p 2 12 /2µ 12 , µ 12 being the reduced mass. The two-fragment model relies on the assumption that the dominant reaction in the quasielastic region consists on the break-up of the nucleus into a knocked-out nucleon with  momentum p fr N and a remaining (A − 1) system in its ground state p fr X , respectively. The relative and center-of-mass momenta of the nucleon and spectator system are obtained as where M X and µ are the mass of the residual nucleus and the reduced mass, respectively. The relative momentum p fr f can be computed in a relativistic way utilizing the correct definition of the final hadronic energy and used to determine the relativistically "fake" kinetic energy (p fr f ) 2 /2µ entering in the energy conserving δ-function of Equation (3). A more detailed discussion of the approach can be found in Rocco et al. [89] and Efros et al. [24].
The results computed within the two-fragment model are displayed by the dashed curves in Figure 3. Note that, now the position and the strength of the LAB rel and ANB rel curves are the same. This implies that the position of the quasielastic peak in the electromagnetic responses no longer depends upon the reference frame and coincides with that of the ANB frame (solid blue curve).
In order to compute the inclusive electromagnetic cross section of Equation (4), R L and R T have to be evaluated for several values of ω and |q|. Hence, due to the sizable computational effort associated with the inversion of the Euclidean response for a given value of |q|, the direct evaluation of Equation (4) is not feasible within GFMC. In order to overcome this limitation, an interpolation procedure based on the concept of scaling was devised in Rocco et al. [24] and used for an accurate and efficient interpolation of the GFMC responses. The right panel of Figure 3 displays the double-differential electron-4 He cross sections for E e = 1108 MeV and θ = 37.5. The red and blue lines have been obtained including the one-body and one-plus two-body contributions in the electromagnetic currents. The inclusion of the meson exchange current leads to a substantial enhancement in the quasielastic region. The same feature is present in the results of Lovato et al. [23,25] where the electromagnetic responses of 4 He and 12 C have been calculated and compared with the experimental data. Separating the longitudinal and transverse channel it has been observed that the two-body term has a predominantly transverse nature; its contribution is almost vanishing in the longitudinal channel when electromagnetic processes are considered. The red solid line indicates one plus two-body current results obtained in the ANB frame, employing the two-body fragment model to account for relativistic kinematics. The inclusion of relativistic effects leads to a shift of the peak position toward lower values of ω and to a reduction of its width.
Recently, GFMC calculations of the neutral-current responses and cross sections for neutrino scattering off 12 C have also been performed [26]. A description of the two-body charge and current operators used in this work is provided in Lovato et al. [23] and Shen et al. [40] and references therein. The left panel of Figure 4 shows the different spatial components of the neutral-current response functions in 12 C at momentum transfer |q| = 570 MeV. The vector and axial contributions are shown separately in all cases but for R xy the entire strength is given by the axial-vector interference. The dashed and solid lines have been obtained including the one-body and one-and two-body currents, respectively. Note that also in this case the two-body term significantly increases in magnitude the response functions in the quasielastic region. At variance with the electromagnetic results, the axial component of the two-body operator in the weak neutral charge produce substantial excess strength in the longitudinal channels as clearly visible in the upper figure. While in the transverse response, shown in the middle panel, we have a two-body enhancement both from the axial and vector terms. The xy response function in the lower panel which arises solely on account of this interference is also modified by the two-body contribution. In the right panel of Figure 4 the CC responses in 4 He at momentum transfer |q| = 600 MeV are displayed. In analogy the NC case, we see that the enhancement given by the inclusion of two-body currents is present in all the different channels. However, in the longitudinal one it is not as significant as for the NC scattering. This has to be ascribed to the relative strength of the axial contribution to the total response which in this case is much smaller than the vector one.

SHORT TIME APPROXIMATION
A novel approach to calculate the short-time propagation resulting from two-nucleon dynamics has been recently developed [28]. The STA method utilizes QMC techniques to evaluate path integrals of one-and two-nucleon currents in real time and predict the response function of nuclei in the quasielastic region. The expression of this response is reported in Equation (3) and can be rewritten as where the sum over final states has been replaced with a realtime propagator. In the following, we drop the Lorentz indices to simplify the notation and replace J(q, ω) → J(q) with ω = ω qe , the one-and two-body current operators utilized are the same as for the GFMC results [40]. The main assumption underlying the STA is that only the active pair of nucleons propagate, this qualitatively amounts to rewrite the final state as where |φ 2 f ′ is the correlated two-nucleon state. To evaluate the response function, two completeness relations on the coordinate states are inserted, yielding In the STA the current-current correlator is rewritten keeping the one-and two-body terms the contributions with three or more active nucleons have been neglected. This amounts to include only two-nucleon interactions in the Hamiltonian; the A-nucleon particle propagator is approximated as where H cm ij = P 2 ij /(4m N ) and H rel ij = p 2 ij /m N + v ij . The A − 2 nucleons are treated as static spectators and their energy is assumed to be peaked around a constant valueĒ A−2 .
We put the two equations together and define j µ † The two-nucleon propagator r 12 |e −iH rel 12 t |r ′ 12 is obtained by summing over the bound and continuum eigenstates of H rel 12 . Note that the interaction effects in the active pair are exactly accounted for. It is convenient to rewrite the STA response in terms of an integral over the relative-and center-of-mass energy of the pair as Within the factorization scheme outlined above, interaction effects at the two-nucleon level are fully retained, and the interference between one-and two-body terms in the electromagnetic current operator are consistently accounted for. The low-energy properties of the system (discrete transition and collective excitation of the nucleus) and the correct energy threshold for the quasielastic region can not be described within the STA. For this reason, some corrections have to be introduced in order to recover the exact value of the threshold for the quasielastic scattering. In particular, the response density is folded with a gaussian kernel where ω(e) = e 2 + ω 2 th exp −e/ω and N is defined requiring that the sum rules are preserved The two parameters controlling the shift and width of the gaussian folding are ω th andω, respectively. The values chosen A comparison between the STA and GFMC electromagnetic response function of 4 He is shown in Figure 5 for |q| = 300 and 500 MeV. The dashed line displays the STA results without any knowledge of the threshold while in the full red line the correct behavior at threshold has been enforced as explained in Equation (53). Including these corrections leads to a shift of the response toward larger values of the energy transfer and a redistribution of the strength; while for |q| = 300 MeV this effect is sizable the results at |q| = 500 are only slightly modified. In both configurations the STA results which include the shift accurately reproduce the GFMC ones.
The role played by final-state interactions within the pair is analyzed in the left panel of Figure 6 where the solid and dashed lines correspond to the transverse electromagnetic response of 4 He obtained with and without interactions effects, respectively. Their inclusion leads to a visible shift in the position of the quasielastic peak toward left. The breakdown of the response into one-body current diagonal and off-diagonal terms, interference between one-and two-body currents, and two-body currents only, is also shown. It is interesting to note that the off-diagonal terms, routinely neglected in the IA scheme, provide a negative contribution depleting the response strength. In analogy with the GFMC findings, the interference between one-and two-body currents provides an important enhancement in the quasielastic region; contrary to the pure two-body current which does not provide a significant contribution to the response function in this kinematics. This is likely to be ascribed to the static limit adopted in both the GFMC and STA approaches to derive the non-relativistic expression of the current. As discussed in Dekker et al. [72] in the static limit all energy dependence of the propagator disappears and the resonance behavior in the dip region is not present.
The isospin dependence of pairs in the back-to-back kinematics, i.e., pairs with low initial center-of-mass momentum and high relative momentum, is studied in the right panel of Figure 6. These pairs can be singled-out in the response densities by requiring the pair center-of-mass momentum P being close to |q| and the relative momentum in the final state being large. Figure 6 displays the response densities at fixed energy E cm ≃ |q| 2 /(4m N ) as a function of the relative energy of the pair e. The solid and dot-dashed curves have been obtained with and without interactions within the pair. The comparison between the full results (black line) and the one-body total (magenta line) shows that the two-nucleon currents do not provide a large contribution at low relative energies. While, the full result becomes substantially larger than the one-body current term for e ≥ 250 MeV. The back-to-back momentum distributions of np pairs are known to dominate over pp or nn pairs at high relative momenta because of tensor correlations as discussed in Schiavilla et al. [90]. In addition to that, the two-nucleon currents are almost entirely in the np pairs, and increase the response by roughly a factor of ∼ 2 around at e = 300 MeV.

EXTENDED FACTORIZATION SCHEME
For large values of ω and |q|, a non-relativistic calculation of the hadron tensor is no longer reliable. Since both the final state and the transition currents depends upon the momentum transfer, relativity has to be properly accounted for. Factorizing the hadronic final state allows one to overcome these difficulties by providing a relativistic description of | f and of the current operator. In addition to that, realistic spectral functions are used to accurately model the dynamics of the target nucleus and account for correlation effects. We start by only retaining the one-body current terms and rewriting the hadronic final state as where |p is a plane wave describing the propagation of the final-state nucleon with momentum |p| and energy e(p) =  system. The energy and momentum of the latter are obtained by energy and momentum conservation relations The incoherent contribution to the one-body hadron tensor can be easily obtained from Equation (55) and by inserting a single-nucleon completeness relation where m N is the rest mass of the initial nucleon. We introduced ω defined asω = ω − E + m N − e(k). To describe the scattering off a bound nucleon with momentum k, the four-momentum transfer employed in the hadronic tensor is replaced by q = (ω, q) →q = (ω, q). The factors m N /e(k) and m N /e(k + q) ensures the implicit covariant normalization of the nucleon quadri-spinors. The hole-spectral function P h (k, E) provides the probability distribution of removing a nucleon with momentum k from the target nucleus, leaving the residual (A − 1)-nucleon system with an excitation energy E. The calculation of the spectral function of finite nuclei is a challenging problem that has seen the endeavor of multiple theory groups. In this work we will focus on the results obtained within the Correlated Basis Function and the Self Consistent Green's Function many-body methods, shortly outlined in section 5.1. For low and moderate values of |q|, interactions between the struck particle and the spectator system become relevant. For this reason, the IA results must be modified to include them [30]. This is achieved by including in the energy spectrum of the propagating nucleon the real part of the optical potential U of Cooper et al. [91] which accounts for its interactions with the mean-field created by the residual system. This potential has to be evaluated for a given kinetic energy of the nucleon t kin (p) = p 2 + m 2 − m, and modifies its energy as e(k + q) = e(k + q) + U t kin (k + q) .
The rescattering processes of the propagating nucleon are described by a convolution scheme which amounts to fold the IA responses with a function f k+q , normalized as The one-body hadron tensor then reads where a generalization of the Glauber theory is utilized to derive the folding function [92] f A more detailed discussion on how to obtain the nuclear transparency T p and the finite width function F p (ω) can be found in Benhar et al. [29] and Benhar et al. [93]. The inclusion of two body-currents requires an extension of the factorization ansatz of Equation (55). In Benhar et al. [34] and Rocco et al. [35,36] the amplitudes involving two-nucleon currents have been included by rewriting the hadronic final state as where |p p ′ a = |p p ′ − |p ′ p is the anti-symmetrized twonucleon plane wave state. The two-body current component of the hadron tensor reads [34] W µν 2b (q, ω) = where P h (k, k ′ , E) is a two-hole spectral function which in Benhar et al. [34] and Rocco et al. [35,36] has been approximated as the product of the one-nucleon ones. Note that, while this is correct for infinite nuclear matter, its application to medium mass nuclei such as 12 C is questionable and should be further investigated.
The production of real pions in the final state will be crucial for the correct understanding of the DUNE results. In order to include this reaction mechanism, we can write the hadronic final state as where p π denotes both the four-momentum (p 0 π , p π ) and the isospin t π of the emitted pion. In analogy with the one-body case reported in Equation (57), the one-body one-pion (1b1π) incoherent contribution to the hadron tensor is given by where e π (p π ) = p 2 + m 2 π is the energy of the outgoing pion. The expression ofω is the same as for the one-body current process. In order to describe the real emission of a pion we need the transition amplitude between the initial one-nucleon state |k to the pion-nucleon |p π p final state. These matrix elements have been obtained within the sophisticated dynamical couple-channel (DCC) model able to describe the πN → πN, γ N → πN, and N(e, e ′ π)N reactions accounting for mesonbaryon channels and nucleon resonances up to an invariant of W = 2 GeV. About 26,000 data points of the πN, γ N → πN, ηN, K , K data from the channel thresholds to W ≤ 2.1 GeV have been fitted to obtain the parameters adopted within the DCC model. The extension to the electroweak sector has been recently performed [38]. Within the DCC model the following Hamiltonian is defined to generate the matrix element p π p|j ν i |k of Equation (65). In the above equation H 0 is the free Hamiltonian while the production of an N * state from a meson-baryon channel c is described by the vertex Ŵ N * ,c . The energy independent mesonexchange potentials v c,c ′ -where c, c ′ = γ N, πN, ηN, K , Kare derived from phenomenological Lagrangians by using the unitary transformation method [94,95]. This hamiltonian is used to generate p π p|j ν i |k of Equation (65). Convoluting the DCC elementary current matrix elements for π-production with the spectral function formalism allows to predict electroweak interactions of finite nuclei in large energy transfer region presented in section 5.2.

Determination of the Hole Spectral Function of Finite Nuclei
The accurate determination of SFs suitable to encompass both single-particle aspects and short-range dynamics is crucial for the theoretical description of lepton-nucleus scattering. Its definition can be given either in terms of nuclear overlaps or as the imaginary part of a two-point hole Green's Function Within the CBF, the hole SF of finite nuclei is written as a sum of two terms [31], P h (k, E) = P MF h (k, E) + P corr h (k, E) displaying distinctly different energy and momentum dependences. The first term is associated to the low momentum and removal-energy region. The spectroscopic factor obtained from (e, e ′ p) scattering measurements are utilized to obtain the first term within a modified mean field (MF) picture. The correlated contribution P corr h (k, E) includes the unbound states of the A − 1 spectator system in which at least one of the spectator nucleons is in a continuum state. The local density approximation (LDA) is adopted to compute this term for finite nuclei where the correlation component of the SF obtained within the CBF theory for isospin-symmetric nuclear matter for a given density ρ is convoluted with the density profile of the nucleus ρ A (R). The applicability of the LDA in obtaining the correlation part of the SF in finite nuclei relies on the observation that short-range nuclear dynamics is not affected by surface and shell effects. This strength of P corr h (k, E) is concentrated in the high momentum and removal energy region as opposed to P MF h (k, E). For momenta larger than the Fermi one, the spectral function coincides with the correlation term. The SCGF approach is a many-body approach which scales polynomially with the number of particle and allows reach nuclei with A up to ∼100. The spectral function is determined within an ab initio theory starting from individual interactions among the nucleons. The central quantity of the SCGF formalism is the one-body Green's function which is directly related to the spectral function through P h (k, E) = − 1 π Im G h (p, p; µ − E) , as expressed in the last equality in Equation (67). This is obtained by adopting an iterative procedure to solve the associated Dyson equation [96,97] where the irreducible self-energy -which encodes nuclear medium effects in the particle [98]-explicitly depends on the propagator itself. To extend the predictive power of the approach to open shell nuclei, the SCGF has been recently reformulated within Gorkov's theory. The particle number is no longer conserved in this new formulation of the propagator in which a grand canonical Hamiltonian is utilized. Breaking of the particle-number symmetry allows one to include pairing correlations and eliminate the degeneracies that would otherwise prevent microscopic calculations for open-shell systems.
The results obtained for the neutron and proton spectral functions of 40 Ar, 40 Ca, and 48 Ti isotopes have been recently presented in Barbieri et al. [99]. The SCGF calculations are performed employing a spherical harmonic oscillator basis in a model space of 14 major shells and varying the frequencȳ h to study the uncertainties resulting from the truncation of the model space. The saturating chiral interactions at next to next to leading order (NNLO sat ) are utilized in order to correctly reproduce radii as well as charge density distributions of nuclei. In the left panel of Figure 7 the charge density profiles computed within the SCGF (solving the Gorkov's equations) or 40 Ca and 40 Ar are compared to experimental data from Emrich et al. [100] and Ottermann et al. [101]. The shaded area in the theoretical curves has been obtained by performing the calculation at the extremes of the rangeh = 14 − 20 MeVwhich from the analysis of Somá et al. [102] turns out to be the optimal one for the convergence of radii and energiesand taking the differences between the results. The authors of Barbieri et al. [99] interpreted this band as a conservative estimate for the theoretical errors due to model space convergence. The right panel of Figure 7 displays the computed hole (particle) P h(p) (p, E) spectral function for neutron removal (addition) from 40 Ar.
In analogy with P MF h (k, E) introduced in the discussion of the CBF results, the peaks present at low E and k correspond to nucleons that occupy the valence shell close to the Fermi surface. In the high momentum and removal energy region, which is typically associated with short range correlation physics, the SCGF spectral function presents a mild tail (not shown in Figure 7). In this regards, it has to be noted that the CBF spectral function relies on the semi-phenomenological AV18 Hamiltonian, which naturally encompass short-range correlations. On the other hand, the NNLO sat interaction is a relatively soft interaction, with a cutoff of 450 MeV which is able to produce tails for large values of the momentum. However, the strength in that region is significantly smaller than the one obtained using AV18 [103].

Results
In this section we present different scattering results obtained using the CBF and SCGF spectral function. In the left panel of Figure 8 we gauge the differences between the two spectral functions by comparing the results obtained for the doubledifferential cross section of electron-12 C scattering at E e = 620 MeV and θ e = 36 • . In the theoretical results we focused on the quasielastic region, including only the onebody current operator of Equation (17). The dashed and solid  curve correspond to the CBF and SCGF SFs, the blue and red lines have been obtained with and without including FSI effects. Calculations carried out employing the two different many-body approaches are in very nice agreement, although they are obtained from different nuclear interactions. FSI effects have been introduced following the procedure discussed in Equation (60). The overall effect is a shift in the position of the quasielastic peak to the left and a redistribution of the strength which leads to a correct reproduction of the experimental data. The right panel of Figure 8 shows the inclusive electron scattering on 40 Ar at the energy and kinematics of the E12-14-012 JLab experiment compared with the SCGF results with and without FSI displayed by the dot-dashed blue and solid red curve, respectively. The real part of the 40 Ca optical potential taken from Cooper et al. [106]-the one of Ar in not available in the literature-and the folding function of Benhar et al. [93] were adopted to obtain the FSI corrections. The prediction based on the NNLO sat interaction and SCGF spectral function slightly underestimates the experimental data at the quasielastic peak. Overall, there is a small discrepancy which is compatible with the errors associated to the nuclear forces [54].
The left panel of Figure 9 shows the double-differential electron-12 C cross sections for E e = 730 MeV, θ e = 37 • . The theoretical results have been obtained using the CBF spectral function and correcting for FSI effects in the quasielastic region. The solid black line corresponds to the total cross section obtained summing up the contributions associated with the different reaction mechanisms. The dashed blue line displays the one-body current contribution while MEC leading to twonucleon emission are given by the short-dashed red line. The dot-dashed magenta line corresponds to the emission of a real pion and a nucleon. A good agreement between theory and data is observed for the different kinematics analyzed. Note that, the strength associated with pion production is necessary to correctly reproduce the peak in the -production region.
The interference between one-and two-body currents has not yet been included. Although the authors of Benhar et al. [34] argue that the inclusion of this contribution within the factorization scheme leads to a small enhancement the dip region, the GFMC and STA calculations presented in sections 3 and 4 display a significant increase in the transverse response due to the interference contribution. Therefore, the consistent implementation of this term within the spectral function formalism is a necessary step to be undertaken in order to properly compare with the GFMC and STA results.
The results obtained for the double-differential CC ν µ -12 C scattering cross sections are shown in the right panel of Figure 9 for E ν = 1 GeV, θ µ = 30 • . The calculations have been carried out within the same framework employed in the electromagnetic case, utilizing the CBF spectral function and including an axial term in all the current operators. In order to compare with experimental data a folding with the energy distribution of a given neutrino flux is needed.

CONCLUSIONS
The success of current-and next-generation of neutrinooscillation experiments strongly depend on the availability of accurate nuclear physics calculations of the dynamics and electroweak interactions of nuclei with quantified theoretical uncertainties. This motivated the advent of many recent theoretical works focused on improving the description of lepton interactions with nuclei. In this review we outlined the main features of three different many-body approaches.
The GFMC is an ab-initio method which provides extremely accurate predictions for the electroweak response functions of nuclei up to 12 C in which correlations are fully retained. In these results the strength and energy-dependence of two-nucleon processes induced by correlation effects and interaction currents provide a sizable contribution in the quasielastic region. The main limitations of this method come from its non-relativistic and fully inclusive nature, i.e., the transition to a given hadronic final state can not be easily identified. Choosing a reference frame that minimizes nucleon momenta allows to extend the applicability of GFMC to larger lepton energies, correctly reproducing experimental data for a large set of kinematics.
For moderate values of momentum transfers, comparing the predictions based on more approximate schemes of nuclear dynamics with the accurate GFMC results in the quasielastic region is extremely important in order to test and validate them.
In the STA, two-nucleon physics is fully accounted for including ground-state correlations and final state interactions among the pair. Within this approach the interference between one-and two-nucleon current is included. When the physical threshold is enforced in the STA calculations, a very good agreement with the GFMC quasi-elastic response of light nuclei is observed for momentum transfers near and above the Fermi momentum.
The formalism based on IA and realistic SF combines a non-relativistic description of the target nucleus including realistic interactions with relativistic currents and kinematics. The results obtained utilizing two-different spectral functions corresponding to the CBF and SCGF calculations have been compared in the quasielastic region. The original formulation of the factorization scheme utilized in the IA only included one-nucleon matrix elements, its recent generalization allowed to include meson-exchange currents and pion-production mechanisms. The elementary amplitudes corresponding to pion-production processes were computed capitalizing on the sophisticated DCC model [37][38][39], which provides robust predictions up to an invariant mass of W ≤ 2.1 GeV. Theoretical calculations that include different reaction mechanisms are in good agreement with inclusive electron-scattering off 12 C.
Over the last decade, we witnessed a great progress in the development of many-body techniques aimed at studying nuclear properties and interactions. These methods rely on nuclear EFT to consistently derive the nuclear Hamiltonian and many-body currents. Despite this success, their application on the broad energy region explored by oscillation experiments involves non-trivial difficulties that have to be explored and understood, such as how to transition to regions where resonance-production occur and how to correctly include relativistic effects. To address these points, models based on a factorization of the hadronic final state have been developed. Benchmark these models with the low-energy description provided by correct many-body calculations and understand the limit of validity of the approximations done will play a crucial role in providing robust predictions of neutrino-nucleus interactions.

AUTHOR CONTRIBUTIONS
The author made a substantial, direct and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
The author would like to thank her collaborators C. Barbieri, O. Benhar, T. S. H. Lee, W. Leidemann, A. Lovato, S. Nakamura, G. Orlandini, S. Pastore, and V. Soma' for their contributions to the studies presented in this work and for many valuable discussions.