Ab initio study of $\boldsymbol{(\nu_\ell,\ell^-)}$ and $\boldsymbol{(\overline{\nu}_\ell,\ell^+)}$ inclusive scattering in $^{12}$C: confronting the MiniBooNE and T2K CCQE data

We carry out an ab initio calculation of the neutrino flux-folded inclusive cross sections, measured on $^{12}$C by the MiniBooNE and T2K collaborations in the charged-current quasielastic (CCQE) regime. The calculation is based on realistic two- and three-nucleon interactions, and on a realistic nuclear electroweak current with one-and two-nucleon terms that are constructed consistently with these interactions and reproduce low-energy electroweak transitions. Numerically exact quantum Monte Carlo methods are utilized to compute the nuclear weak response functions, by fully retaining many-body correlations in the initial and final states and interference effects between one- and two-body current contributions. We employ a nucleon axial form factor of the dipole form with $\Lambda_A = 1.0$ or $1.15$ GeV, the latter more in line with a very recent lattice QCD determination. The calculated cross sections are found to be in good agreement with the neutrino data of MiniBooNE and T2K, and antineutrino MiniBooNE data, yielding a consistent picture of nuclei and their electroweak properties across a wide regime of energy and momenta.


I. INTRODUCTION
There is a large program of accelerator neutrino experiments in operation or in the planning phase in the US and elsewhere to measure the parameters that characterize the probabilities for flavor oscillations of these particles-mass differences, mixing angles, and the charge-conjugation and parity violating phase. These experiments do not directly measure oscillation probabilities, of course, but rather event-rate distributions as function of the observed energy E in the detector, schematically where φ α (E ν ) is the flux of neutrinos of flavor α (ν α 's) at the source as function of energy E ν , P (ν α → ν β ; E ν ) is the probability for oscillation of a ν α into a flavor ν β , and σ β (E ν , E) is the ν β -nucleus cross section. The neutrino energy is reconstructed from the tracks in the detector of the outgoing lepton in an inclusive scattering setting, and, additionally, the tracks of final hadrons in a semi-inclusive one. As a consequence, the determination of oscillation parameters depends strongly on neutrino interaction physics, since the interactions observed in the detector result from the folding of the energydependent neutrino flux, energy-dependent cross section, and energy-dependent nuclear (strong-and electroweakinteraction) effects. The appreciation of these difficulties has led, in the last decade or so, to a flurry of activity by nuclear theorists, who have attempted to provide accurate estimates for neutrino-nucleus (ν-A) inclusive (and semi-inclusive) cross sections (for a summary of efforts in this area see Ref. [1]). This is a very challenging task, primarily, because neutrino fluxes in current (such as MiniBooNE, T2K, MicroBooNE, and Minerνa) and future (DUNE) experiments extend over a rather wide energy range, from threshold to, in several cases, multi-GeV energies. Thus, observed ν-A cross sections, resulting from the folding in Eq. (1), may include contributions from energy-and momentum-transfer regions of the nuclear weak response where drastically different dynamical regimes are at play, from the structure and collective behavior of low-lying nuclear excitations in the threshold region, to the quark substructure of individual nucleons in the deep inelastic region. Moreover, for some of the nuclear targets employed in the detectors of these experiments, such as 40 Ar (MicroBooNE), and 56 Fe and 208 Pb (Minerνa), the full structure of the ground states is difficult to calculate exactly.
Theoretical studies have attempted to provide a description of the nuclear weak response in this wide range of energy and momentum transfers. They typically rely on a relativistic Fermi gas (RFG) [2][3][4] or relativistic mean field (RMF) [5][6][7][8] picture of the nucleus. Some, notably those of Refs. [9][10][11][12][13], include correlation effects in the random-phase approximation (RPA) induced by effective particle-hole interactions in the N -N , N -∆, ∆-N , and ∆-∆ sectors, use various inputs from pion-nucleus phenomenology, and lead to predictions for electromagnetic and strong spin-isospin response functions of nuclei, as measured, respectively, in inclusive electron scattering and in pion and charge-exchange reactions, in reasonable agreement with data. Some utilize the phenomenologi-cal SuperScaling (SuSA) approach-scaling with nuclear mass number [14]-in which a universal scaling function, derived from analyses of (e, e ) data on a number of nuclei, is used to obtain estimates for the corresponding ν-A cross sections [15,16]. Recently, SuSA, which has proven to be quite successful, has been extended (SuSAv2) by incorporating elements from the RMF approach to account for differences between the vector and axial components of the weak current, and between their isoscalar and isovector content [17][18][19]. Yet others rely on factorization of the hadronic final state and realistic spectral functions S(p m , E m ) to combine an accurate description of the nuclear ground state with relativistic currents and kinematics. Spectral functions of atomic nuclei are calculated either with microscopic methods-for example, the self-consistent Green's function technique [20][21][22][23]-or by combining inputs from (e, e p) data to characterize the low missing-momentum and missing-energy region with accurate many-body calculations of the nuclear matter spectral function [24][25][26][27][28] folded with the single-nucleon density to describe the "correlation region", corresponding to high missing energies and momenta [29,30].
Many of the above models have achieved a remarkable phenomenological success and much improved agreement with experimental data, compared to simple RFG calculations. However, it is fair to note that they rely on a somewhat approximate description of nuclear dynamics that does not fully capture correlation effects in both the initial and final states and does not generally use as inputs realistic nuclear interactions and consistent electroweak currents. Hence, it is important to carefully assess their validity-especially in the axial sector-by testing them against more microscopic calculations.
In the present study we report on an ab initio calculation of the nuclear weak response induced by chargedcurrent (CC) (ν , − ) and (ν , + ) processes. The strong interactions among nucleons are represented by two-and three-body terms, while their coupling to the electroweak field is accounted for by one-and two-body currents (see the reviews [31,32] and references therein). The twobody interaction [33] is constrained by fits to the nucleonnucleon (N N ) database up to lab energies of 350 MeV (albeit it provides a good description of the N N cross section well beyond the pion production threshold, up to 500 MeV or so). The three body interaction [34] is calibrated by a fit to the energies of a number of low-lying nuclear states in the mass range A = 3-10.
The one-body currents follow from a non-relativistic expansion of the covariant single-nucleon CC, including nucleon electroweak form factors consistent with available experimental data. In particular, results reported in Sec. III are obtained by using a dipole axial form factor with cutoff Λ A equal to either 1.0 GeV or 1.15 GeV. The former has been extracted from proton and deuteron experiments [35][36][37][38], while the latter is obtained by recent lattice QCD calculations [39,40] that also reproduce the vector form factors measured in electron scattering.
Two-body currents are derived from meson-exchange phenomenology including pion and ρ-meson exchanges as well as N -to-∆ transition currents (with ∆'s taken, however, in the static limit) [41]. The short-range behavior of these currents is prescribed to be consistent with that of the two-nucleon interaction [42,43]. In the vector sector, they contain no free parameters, while in the axial sector the single unknown parameter present-the N -to-∆ axial coupling constant-is fixed by reproducing the experimental value of the tritium Gamow-Teller matrix element. The theoretical framework outlined above (and discussed more expansively in Sec. II below) has been shown to provide, in numerically accurate quantum Monte Carlo (QMC) calculations, a quantitatively successful description of a large body of experimental data on light nuclei (A ≤ 12), including, among others, energy spectra of low-lying states, static properties (magnetic and quadrupole moments), low-energy radiative and weak transition rates, electromagnetic ground and transition form factors, and electroweak dynamic response (for a review, see [32] and references therein). Especially relevant in the present context are the QMC studies of the 12 C electromagnetic ground-state structure [44], and longitudinal and transverse response functions at intermediate momentum transfers q in the (300-700) MeV range, and for energy transfers ω in the quasielastic region [45].
However, it is also important to recognize the limitations inherent to the approach we have adopted here: firstly, it addresses only inclusive scattering; secondly, it does not account for explicit pion production mechanisms and therefore cannot describe the nuclear electroweak response in the ∆ resonance region and beyond; and thirdly, it relies on what is in essence a non-relativistic formulation of the dynamics and electroweak currents. 1 These limitations notwithstanding, it should be emphasized that in the quasielastic regime specified earlier, this approach includes all of the relevant physics for inclusive scattering and is expected to be quite accurate. It is for this reason that we compare our predictions (in Sec. III) for the 12 C flux-averaged inclusive cross sections-differential in the outgoing lepton energy and scattering angle-to the MiniBooNE and T2K CC "quasielastic" (CCQE) data sets [46][47][48]. These data sets are characterized by the absence of pions in the final state. Clearly, their interpretation as purely "quasielastic" is complicated by pions that are created at the interaction vertex and are subsequently reabsorbed in the nuclear medium [49]. The unambiguous identification of these contributions is problematic, and modeldependent at best, requiring an accurate modeling of both the pion-production cross section and subsequent reabsorption (and their interference). Currently, they are estimated using Monte Carlo event generators. As a consequence, experimentally-extracted CCQE cross sections are accompanied by significant uncertainties.

II. CALCULATION
The inclusive double-differential cross section for a charged-current scattering process initiated by a neutrino off a nuclear target can be expressed as where the − or + sign corresponds to a neutrino (ν) or antineutrino (ν) induced reaction. We adopt the values G F = 1.1803 × 10 −5 GeV −2 , corrected for the bulk of the inner radiative corrections [50], and cos θ c = 0.97425 [51]. The initial ν (or ν) and final lepton four-momenta are, respectively, k ν = (E ν , k ν ) and k = (E , k ), T is the kinetic energy of the lepton (rest mass m ), and θ is its scattering angle relative to the incoming neutrino direction. The kinematical factors v αβ associated with the contraction of the leptonic tensor, in the general case in which the dependence on m is kept, are reported in Appendix A of Ref. [41]. The nuclear response functions encode all information on nuclear structure and dynamics, and are defined, in a schematic notation, as (see Ref. [41] for explicit expressions) where |i represents the Z A ground state of energy E i , |f represents the bound or scattering state of the final Z+1 A or Z−1 A nuclear system, depending on whether the (ν , − ) or (ν , + ) process is being considered, of energy E f , j α CC (q, ω) are the relevant components of the weak charged current (CC), and an average over the initial spin projections of Z A is understood (note, however, that the 12 C ground state has spin-parity assignments J π = 0 + ). The dynamical framework adopted in the calculations below has been described elsewhere in considerable detail, most recently in the review [32]. Next, we provide a brief description for completeness.

A. Interactions and currents
Strong interactions are described by two-and threenucleon terms, respectively, the Argonne v 18 [33] (AV18) and Illinois-7 [34] (IL7) models. The AV18 reproduced the nucleon-nucleon database available at the time (1995) with a χ 2 /datum close to one [33] for lab kinetic energy up to 350 MeV, slightly above the pion production threshold. Even today that the database has increased in size considerably (to over 5,200 data points over the energy range 0-300 MeV), the AV18 still gives (without a refit) a very respectable χ 2 /datum of about 1.5 [52]. The IL7 three-nucleon interaction model contains a small number (4) of parameters, which characterize the overall strengths of two-and multi-pion exchange terms involving ∆-isobar excitations, and of a purely phenomenological (isospin-dependent) central term. These parameters are constrained by a fit to the energies of about 23 lowlying nuclear states with mass number A in the range 3-10 [53]. The resulting AV18+IL7 Hamiltonian then leads, in accurate QMC calculations, to predictions for about 100 ground-and excited-state energies up to A = 12, including the 12 C ground-and Hoyle-state energies, in good agreement with the corresponding empirical values [32].
Electroweak probes couple to single nucleons (impulse approximation) as well as to clusters of nucleons via oneand many-body currents. The CC model adopted in the present study, identical to that of Ref. [41] and most recently employed to compute the muon-capture inclusive rates on 3 H and 4 He [54], contains one-and two-body terms. The former are derived from the covariant singlenucleon CC in a non-relativistic expansion that retains corrections proportional up to the inverse square of the nucleon mass. Two-body (vector and axial) terms arise from effective π-and ρ-meson exchanges, and N -to-∆ excitations, treated in the static limit. A ρπ transition mechanism is also included in the axial component. In GFMC calculations we utilize configuration-space representations of these currents, regularized by a prescription which, by construction, makes their short-range behavior consistent with the AV18 interaction [31]. The value for the transition (axial) coupling constant g * A in the Nto-∆ axial current is determined by reproducing, within the present dynamical framework, the measured Gamow-Teller matrix element contributing to tritium β-decay, and is listed in Table I (Set I) of Ref. [41], where explicit expressions for these currents can also be found.
The (isovector) nucleon form factors in the CC vector component are taken as functions of the squared fourmomentum transfer (Q 2 = q 2 − ω 2 ) from a modern fit to the available electron scattering data [55] (in contrast to Ref. [41], in which we adopted a simple dipole parametrization of these form factors). The axial form factor G A (Q 2 ) of the nucleon is of a dipole form with a cutoff mass of either 1 GeV or 1.15 GeV, while its induced pseudoscalar form factor, derived from the PCAC constraint and pion-pole dominance, is in accord with values extracted from precise measurements of the muoncapture rate on hydrogen and 3 He [56] as well as with predictions based on chiral perturbation theory [57,58]. Lastly, the N -to-∆ transition form factor in the vector sector is as obtained in an analysis of γN data in the ∆-resonance region [59], while that in the axial sector, because of the lack of available experimental data, is simply taken to have the same functional form of Values for the parameters entering these axial form factors are specified in Ref. [41].

B. Electroweak response functions
The calculation of the response functions in Eq. (3) proceeds in two steps. The first consists in Laplacetransforming R αβ (q, ω) with respect to ω, which reduces to the following current-current correlator (Euclidean response function) where H is the Hamiltonian (here the AV18-IL7 model).
The energy dependence of j α CC (q, ω) comes in via the nucleon and N -to-∆ transition form factors, which are taken as functions of Q 2 , as noted above. We freeze the ω-dependence by fixing Q 2 at the value Q 2 qe = q 2 − ω 2 qe with the quasielastic energy transfer ω qe given by ω qe = q 2 + m 2 − m (m is the nucleon mass). This is needed in order to exploit the completeness over the nuclear final states in evaluating the Laplace transforms of R αβ (q, ω). The correlator is then computed with Green's function Monte Carlo (GFMC) methods [44,45,[60][61][62][63].
It should be stressed that no additional approximations are made beyond those inherent to the modeling of the interactions and currents. The response is thus calculated ab initio by treating completely correlations in the initial state, by accounting consistently through the imaginarytime propagation for interaction effects in the final states, and, in particular, by retaining in full the important interference between one-and two-nucleon currents.
Because of the computational cost of the present study (of the order of 130 million core hours on the massively parallel computer MIRA at ANL), however, we only propagate the Z−1 A system, i.e., j α CC in Eq. (4) is the charge lowering current corresponding to the process (ν , + ). If electromagnetic interactions and isospinsymmetry-breaking terms in the strong interactions were to be ignored, the final states |f ; Z+1 A and |f ; Z−1 A of the Z+1 A and Z−1 A nuclear systems would simply be related to each other via |f where τ i,x is the isospin flip operator converting proton i into a neutron or viceversa. Matrix elements of the charge-raising and charge-lowering current between the Z A state and, respectively, the Z+1 A and Z−1 A states would then be identical. We will assume here this is the case for 12 C, and obtain the response functions corresponding to the (ν , − ) process from those corresponding to the (ν , + ) process by correcting the final state energies of the 12 B system by the difference in ground-state energies between 12 N and 12 B-in practice, by shifting the response functions by about 5.5 MeV. We expect this approximation to be inaccurate in the threshold region; however, in quasielastic kinematics and beyond, it should be of little import. The second step employs maximum-entropy techniques, developed specifically for this type of problem in Ref. [62] (a fairly complete account of them is given in that work), to perform the analytic continuation of the Euclidean response functions, corresponding to the "inversion" of the Laplace transforms. The resulting R αβ (q, ω qe ) are rescaled as follows to account for the correct ω-dependence of the various form factors. The 00, 0z, zz, and xx response functions are given by the incoherent sum of the (squared) matrix elements associated with the CC vector (V ) and axial (A) components, while the xy response function involves interference between these components. The 00, 0z, and zz V contributions are multiplied by the factor G V and G V M are the isoscalar and isovector combinations of the proton (p) and neutron (n) electric (E) and magnetic (M ) form factors (in the parametrization of Ref. [55]). These multiplicative factors naturally emerge by considering the dominant one-body terms in the CC V current.
The 0, 0z, zz, and xx A contributions are multiplied by the factor For these contributions such a rescaling turns out to fully restore the correct ωdependence, since the one-and two-body axial currents, including those associated with ∆-isobar intermediate states, are proportional to G A (Q 2 ) in the present modeling [41]. Lastly, the interference response is rescaled by the factor . Below, we show that the procedure above essentially accounts for the correct ω-dependence implicit in the complete CC response.
The five response functions entering the CC cross section have been calculated with GFMC methods for momentum transfers in the range (100-700) MeV in steps of 100 MeV. To reduce clutter, we present in Fig. 1 only those obtained at q = 300, 500, and 700 MeV (note that the scales for R αβ are different in each panel). 2 The transverse (xx) and interference (xy) response functions are largest but of opposite sign (the xy response as defined here is negative). Consequently, the contributions v xx R xx and v xy R xy in the CC cross section add up for neutrino scattering and tend to cancel each other out for antineutrino scattering (the kinematical factors v xx and v xy are positive [41]).
Two-body terms in the CC significantly increase the magnitude of the response functions obtained in impulse approximation (i.e., with one-body currents), over the whole quasielastic region, except for R 00 at low ω. This increase in strength mostly comes about because of constructive interference between the one-and two-body current matrix elements, and is consistent with that expected on the basis of sum rule analyses [64]. Two-2 Tabulations of GFMC-calculated R αβ (q, ω) for q in the range (100-700) MeV and ω from threshold to ω q are available upon request. body contributions are found to be especially largeaccounting for more than 50% of the total calculated strength-in R zz , which involves the longitudinal components (along the direction of the three-momentum transfer) of the CC.

C. Scaling analysis
The analysis of scaling properties of nuclear response functions has proven to be a useful tool to elucidate important aspects of the many-body dynamics in the quasielastic region. Scaling occurs when the electroweak response functions, divided by appropriate pre-factors describing single-nucleon physics, no longer depend upon the momentum q and energy transfer ω, but only on a specific function of them ψ(q, ω), yielding where k F is the Fermi momentum of the system. In the non-relativistic limit, the scaling variable is given by [65] where is introduced to account for nuclear binding effects.
The pre-factors associated with the electromagnetic longitudinal and transverse responses can be found in Ref. [65]. Here we extend the scaling analysis to the five response functions relevant for neutrino-nucleus scattering induced by CC transitions. The (longitudinal and transverse) pre-factors associated with vector currents are related to those of (isovector) electromagnetic currents by the CVC constraint; the pre-factors associated with axial currents bring about additional terms, whose relativistic expressions can be found in Ref. [66].
Within the Fermi gas model [3], the following scaling function can be analytically derived, by assuming one-body currents only. However, unlike the latter expression, which is symmetric and centered around ψ = 0, the scaling functions extracted from experimental data and those inferred from more realistic models of nuclear dynamics exhibit a clearly asymmetric shape, with a tail extending in the ψ > 0 region [67]. Moreover, while the Fermi gas scaling function is universal and does not depend upon the specific transition operator, such is not the case when the spin and charge dependence of nuclear interactions in the final states are taken into account [61]. The xx (xy) scaling functions displayed in the upper two (lower two) panels of Fig. 2 have been obtained as in Eq. (5), i.e., by dividing the GFMC electroweak response functions in the transverse (interference) channel by the appropriate pre-factors for the CC vector and axial components. The upper and lower panels for each set (xx and xy) correspond to including one-body only, and one-and two-body, current operators. The dotted (red), dashed (blue), and solid (green) lines show the xx and xy scaling functions for q = 500, 600, and 700 MeV, respectively. The shaded area indicates the uncertainty in the maximum-entropy inversion procedure and also reflects the statistical errors of the GFMC calculations. The xy scaling functions, shown in the lower two panels of Fig. 2, are almost identical to the xx ones in both cases (one-body only, and one-and two-body currents).
In contrast to the Fermi gas model, nuclear correlations in the initial and final states, which are exactly treated in the GFMC method, yield asymmetric scaling functions, with tails that extend well beyond ψ > 1. Note that the scaling functions can be significantly different in the other channels, for example the longitudinal and transverse response in electron scattering [27].
The different curves clearly exhibit a scaling behavior, as they are almost independent of momentum transfer. This is expected to be even more accurate at larger q values. More interesting is the observation that scaling persists even when two-body current contributions are included in the response functions, as shown in the second and fourth panels of Fig. 2. While these contributions generate significant excess strength in f xx and f xy , they do not spoil their scaling properties. An explanation of these features can be found in Ref. [68] for the case of the electromagnetic response, and similar considerations remain valid here. In essence, in the xx and xy responses the excess strength seen in the quasielastic region comes about because two-body currents lead to final states which are very similar to those produced by an electroweak interaction vertex on a single nucleon followed by the subsequent high-momentum strong interaction of this nucleon with another nucleon. The resulting (constructive) interference between the corresponding matrix elements generates excess strength which is spread out over the quasielastic peak region in a way very similar to the response arising from the high-momentum part of the single-nucleon currents associated with pion exchange interactions. We defer to Ref. [68] for a more comprehensive discussion of scaling in the present context of microscopic Hamiltonians and currents. This reference also discusses superscaling [14]-scaling with respect to the mass number-and the absence of scaling observed in the ∆-resonance region.
An analogous scaling behavior is also seen in the 00, 0z, and zz channels. We capitalize on this feature in order to extrapolate the response functions at large momentum transfers q > q = 700 MeV, that is, beyond the range of those calculated with GFMC methods. It turns out they are needed when computing flux-folded cross sections (see Sec. III below). We parametrize them as where f αβ (ψ) are the scaling functions determined from the GFMC-calculated responses at q. The underlying assumption is that the f αβ (ψ) for q > q coincide with those at q. To account for the small scaling violations, we conservatively associate an uncertainty to this extrapolation procedure corresponding to twice the difference between the scaling functions at q = 600 and 700 MeV.

III. RESULTS
Muon neutrino and antineutrino flux-averaged cross sections are obtained from where φ(E ν ) is the normalized ν µ or ν µ flux-those for MiniBooNE and T2K are shown in Fig. 3-and dσ(E ν )/(dT µ d cos θ µ ) are the corresponding inclusive cross sections of Eq. (2). The experimental data are binned in cos θ µ bins of constant width (0.1) for Mini-BooNE, and varying widths for T2K; when comparing to these data, the calculated cross sections are averaged over the relevant cos θ µ bin. Predictions for the flux-averaged cross sections on 12 C corresponding to the two experiments and obtained by including one-body only, and one-and two-body, currents are shown by, respectively, dashed (green) and solid (blue) lines in Figs. 4-6. The shaded areas result from combining statistical errors associated with the GFMC evaluation of the Euclidean response functions, uncertainties in the maximum-entropy inversion of them, and uncertainties due to extrapolation of the response functions outside the calculated (q, ω) range, which is 100 MeV ≤ q ≤ 700 MeV and ω from threshold to ω q. This extrapolation is carried out by exploiting the scaling property of the various response functions, as outlined at the end of the previous section. The large cancellation between the dominant terms proportional to v xx R xx and v xy R xy in antineutrino cross sections leads to somewhat broader error bands than for the neutrino cross sections, for which those terms add up. Furthermore, we note that the cross-section scales in Figs. 4 and 5 are different, those for the ν µ -CCQE data being a factor of about 2 to 10 smaller than for the ν-CCQE data as the muon scattering angle increases from 0 • to 90 • .
Overall, the MiniBooNE ν µ and ν µ , and T2K ν µ , data are in good agreement with theory, when including the contributions of two-body currents. This is especially noticeable in the case of the MiniBooNE ν µ data at forward scattering angles. However, the calculated cross sections underestimate somewhat the MiniBooNE ν µ data at progressively larger muon kinetic energy T µ and backward scattering angles θ µ , and the ν µ data at forward θ µ over the whole T µ range. By contrast, the full theory (with one-and two-body currents) appears to provide a good description of the T2K ν µ data over the whole measured region.
For a given initial neutrino energy E ν , the calculated cross section is largest at the muon energy T µ corresponding to that of the quasielastic peak, where m is nucleon mass, and on the r.h.s. of the equation above we have neglected the muon mass. The position of the quasielastic peak then moves to the left, towards lower and lower T qe µ , as θ µ changes from the forward to the backward hemisphere. The general trend expected on the basis of this simple picture is reflected in the calculation and data, even though the cross sections in Figs. 4-6 result from a folding with the neutrino flux, which is far from being monochromatic. Nevertheless, the correlation between peak location in the flux-averaged cross sections and θ µ remains. For example, the T2K flux is largest at E ν ≈ 560 MeV and fairly narrow; hence, one would expect the T2K flux-averaged cross section be peaked at the muon momentum p qe µ ≈ 550 MeV for cos θ µ = 1, and p qe µ ≈ 450 MeV for cos θ µ = 0.65, in reasonable accord with the data of Fig. 6.
In Figs. 4 and 5 we also present the flux-folded ν µ and ν µ cross sections obtained in plane-wave-impulseapproximation (PWIA) for three different bins in cos θ µ (corresponding to the forward, intermediate, and backward region) of the MiniBooNE data. We have adopted here the most naive (non-relativistic) formulation of PWIA based on the single-nucleon momentum distribution rather than the spectral function. 3 Hence, the PWIA response functions follow from where the factors x αβ (p, q, ω) denote appropriate combinations of the CC components (the same single-nucleon CC utilized in the GFMC calculations), and N (p) is the nucleon momentum distribution in 12 C (as calculated in Ref. [69]). The effects of nuclear interactions are subsumed in the single parameter E, which can be interpreted as an average separation energy (we take the value A couple of comments are in order. First, the cross sections in PWIA are to be compared to those obtained with the GFMC method by including only one-body currents (curves labeled GFMC 1b): they are found to be systematically larger than the GFMC predictions, particularly at forward angles. Furthermore, it appears that the (spurious) excess strength in the PWIA cross sections (in the same forward-angle kinematics) matches the in-crease produced by two-body currents in the GFMC calculations (difference between the GFMC 1b and GFMC 12b curves). This should be viewed as accidental.
Second, the PWIA and PWIA-R cross sections are very close to each other, except in the ν case at backward angles. In this kinematical regime there are large cancelations between the dominant terms proportional to the transverse and interference response functions; indeed, as θ µ changes from 0 • to about 90 • , the ν cross section drops by an order of magnitude. As already noted, these cancellations are also observed in the complete (GFMC 12b) calculation, and lead to the rather broad uncertainty bands in Fig. 5. Aside from this qualification, however, the closeness between the PWIA and PWIA-R results provides corroboration for the validity of the rescaling procedure of the electroweak form factors, needed to carry out the GFMC computation of the Euclidean response functions.  Fig. 4 but for νµ-CCQE scattering. The experimental data and their shape uncertainties are from Ref. [47]. The additional 17.4% normalization uncertainty is not shown here.

IV. CONCLUSIONS
We have reported on an ab initio study, based on realistic nuclear interactions and electroweak currents, of neutrino (and antineutrino) inclusive scattering on 12 C in the CCQE regime of the MiniBooNE and T2K data. Nuclear response functions have been calculated with QMC methods and, therefore, within the description of nuclear dynamics that we have adopted here, fully include the effects of many-body correlations induced by the interactions in the initial and final states, and correctly account for the important (constructive) interference between one-and two-body current contributions. This interference leads to a significant increase in the crosssection results obtained in impulse approximation, and is important for bringing theory into much better agreement with experiment.
The nucleon and nucleon-to-∆ electroweak form factors entering the currents have been taken from modern parameterizations of elastic electron scattering data on the nucleon and deuteron, and neutrino scattering data on the proton and deuteron. In particular, the Q 2dependence of the nucleon axial form factor G A (Q 2 ) is of a dipole form with a cutoff Λ A ≈ 1 GeV. The nucleonto-∆ axial coupling constant g * A has been fixed by reproducing the Gamow-Teller matrix element measured in tritium β decay, while the Q 2 -dependence of its (transition) form factor G * A (Q 2 ) has simply been assumed to be the same as that of G A (Q 2 ), since no experimental information is currently available on G * A (Q 2 ). First-principles LQCD calculations of nucleon (and, possibly, nucleon-to-∆) electroweak form factors could potentially have a significant impact on calculations of neutrino-nucleus cross sections, since these form factors constitute essential inputs to the nuclear CC. This is especially the case for G A (Q 2 ) and the induced pseudoscalar form factor G P (Q 2 ), whose Q 2 -dependence is experimentally poorly known. In this context, it is interesting to note that recent LQCD studies [39,40,70] find the Q 2 fall-off of G A (Q 2 ) with increasing Q 2 significantly less drastic than implied by the dipole behavior with Λ A ≈ 1 GeV. They also find the nucleon isovector vector form factors in agreement with experimental data which are of course quite accurate. These calculations suggest a larger value of Λ A may be appropriate. We investigate the implications of this finding by presenting in Fig. 7 the flux-folded cross sections (for MiniBooNE and selected bins in cos θ µ ), obtained by replacing in the dipole parametrization the cutoff Λ A ≈ 1 GeV with the value Λ A ≈ 1.15 GeV. As expected, this leads generally to an increase of the GFMC predictions over the whole kinematical range. Since the dominant terms in the cross section proportional to the transverse and interference response functions tend to cancel for ν µ , the magnitude of the increase turns out to be more pronounced for ν µ than for ν µ -as a matter of fact, the ν µ cross sections are reduced at backward angles (0.1 ≤ cos θ µ ≤ 0.2). Overall, it appears that the harder cutoff implied by the LQCD calculation of G A (Q 2 ) improves the accord of theory with experiment, marginally for ν µ and more substantially for ν µ . In view of the large errors and large normalization uncertainties of the MiniBooNE and T2K data, however, we caution the reader from drawing too definite conclusions from the present analysis. Indeed more precise nucleon form factors can be obtained through further lattice QCD calculations or experiments on the nucleon and deuteron, respectively.
Of course, many challenges remain ahead, to mention just three: the inclusion of relativity and pion-production mechanisms, and the treatment of heavier nuclei (notably 40 Ar). While some of these issues, for example the implementation of relativistic dynamics via a relativistic Hamiltonian along the lines of Ref. [71], could conceivably be incorporated in the present GFMC approach, it is out of the question that such an approach could be utilized to describe the ∆-resonance region of the cross section or, even more remotely, extended to nuclei with mass number much larger than 12, at least for the foreseeable future. In fact, it maybe unnecessary, as more approximate methods exist to deal effectively with some of these challenges, including factorization approaches based on one-and two-nucleon spectral functions [28,72] or on  [39]. The first two rows correspond to the MiniBooNE flux-folded νµ and νµ CCQE cross sections, respectively; the last row corresponds to the T2K νµ CCQE data. In the theoretical curves the total one-plus two-body current contribution to the cross section is displayed. the short-time approximation of the nuclear many-body propagator [68] for relativity and pion production, and auxiliary-field-diffusion Monte Carlo methods [73] to describe the ground states of medium-weight nuclei. We are optimistic that the next few years will witness substantive progress in the further development and implementation of these approximate methods to address the high-energy region of the nuclear electroweak response.
Finally, factorization approaches can also be helpful in obtaining some information on exclusive final states. For more complete treatment of these or, in fact, low energy peaks in the threshold region of the response quantum computers could play a role, given sufficient size and sufficiently low error rates [74,75].