Lattice QCD and Neutrino-Nucleus Scattering

TL;DR

This paper explores how lattice QCD calculations of nucleon and nuclear structures can improve the understanding of neutrino-nucleus interactions, crucial for current and future neutrino experiments like DUNE and HyperK.

Contribution

It discusses the potential and importance of lattice QCD in providing nuclear matrix elements for neutrino physics, linking theory with experimental needs.

Findings

Lattice QCD can compute nucleon and nuclear matrix elements relevant to neutrino interactions.

These calculations are vital for interpreting data from upcoming neutrino experiments.

The paper highlights the synergy between lattice QCD and nuclear models in neutrino physics.

Abstract

This document is one of a series of whitepapers from the USQCD collaboration. Here, we discuss opportunities for lattice QCD in neutrino-oscillation physics, which inevitably entails nucleon and nuclear structure. In addition to discussing pertinent lattice-QCD calculations of nucleon and nuclear matrix elements, the interplay with models of nuclei is discussed. This program of lattice- QCD calculations is relevant to current and upcoming neutrino experiments, becoming increasingly important on the timescale of LBNF/DUNE and HyperK.

Tables

Table 1. Table 1: Sample of calculations of nucleon form factors going on worldwide. In the first column, “2”, “2+1”, and “2+1+1” all denote two equal-mass quarks for up and down; the latter two include strange and charm, respectively. The last column indicates work in which USQCD members participate.

Sea quarks	Valence quarks	$N_{\text{ens}}$	$a$ (fm)	$M_{\pi }$ (MeV)	Collaboration	Ref.	USQCD
2 Wilson-clover	same as sea	11	0.06–0.08	150–490	RQCD	Bali et al. (2015)
2 TM clover	same as sea	11	0.09	130	ETM	Alexandrou et al. (2017a)
2 Wilson-clover	same as sea	11	0.05–0.08	190–470	Mainz (CLS)	Capitani et al. (2019)
2+1 overlap	same as sea	14	0.11	290–540	JLQCD	Yamanaka et al. (2018)
2+1 domain wall Blum et al. (2016)	overlap	13	0.08–0.15	170–340	$\chi$QCD	Liang et al. (2018a)	✓
2+1 Wilson-clover	same as sea	11	0.085	146, 135	PACS	Shintani et al. (2019)
2+1 Wilson-clover	same as sea	11	0.05–0.09	200–350	Mainz (CLS)	Ottnad et al. (2018)
2+1+1 HISQ Bazavov et al. (2013)	Wilson-clover	18	0.06–0.12	135–210	PNDME	Gupta et al. (2017)	✓
2+1+1 HISQ Bazavov et al. (2013)	domain wall	16	0.09–0.15	130–400	CalLat	Chang et al. (2018a)	✓
2+1+1 TM clover	same as sea	13	0.09–0.15	140	ETM	Constantinou et al. (2018)	✓
2+1+1 HISQ	same as sea	13	0.09–0.15	135	Fermilab/MILC	Meyer et al. (2016b)	✓

Equations

$${\nu }_{\ell }A\to {\ell }^{-}X,{\bar{\nu }}_{\ell }A\to {\ell }^{+}X,$$

$${\nu }_{\ell }n\to l^{-}p,{\bar{\nu }}_{\ell }p\to l^{+}n,$$

$$⟨f|J_{\nu }|i⟩,⟨f|J_{\mu }^{†}{J}_{\nu }|i⟩,⟨f|\mathcal{O}|i⟩,$$

$$⟨p(p^{\prime })|J_{\mu }^{+}|n(p)⟩$$

$$+{\gamma }_{\mu }{\gamma }_5{F}_A^{\text{CC}}(q^2)+{\gamma }_5\frac{q_{\mu }}{M_N}{F}_P^{\text{CC}}(q^2)+{\gamma }_5\frac{(p^{\prime }+p{)}_{\mu }}{M_N}{F}_T^{\text{CC}}(q^2)]u^{(n)}(p),$$

$$r_E^2\equiv 6{\frac{d{G}_E}{d{q}^2}|}_{q^2=0},r_M^2\equiv \frac{6}{G_M(0)}{\frac{d{G}_M}{d{q}^2}|}_{q^2=0},r_i^2\equiv \frac{6}{F_i(0)}{\frac{d{F}_i}{d{q}^2}|}_{q^2=0},$$

$$F_A(q^2)=\frac{g_A}{(1-q^2/m_A^2{)}^2},$$

$$z(t)=\frac{\sqrt{t_{\text{cut}}-t}-\sqrt{t_{\text{cut}}-t_0}}{\sqrt{t_{\text{cut}}-t}+\sqrt{t_{\text{cut}}-t_0}},$$

$$F(z)=\sum _k{a}_k{z}^k$$

$$F_i^{\text{NC}}$$

$$i$$

$$F_A^{\text{NC}}=\frac{1}{2}(-F_A^{\text{CC}}+F_A^s).$$

Full text

EXECUTIVE SUMMARY

In 2018, the USQCD collaboration’s Executive Committee organized several subcommittees to recognize future opportunities and formulate possible goals for lattice field theory calculations in several physics areas. The conclusions of these studies, along with community input, are presented in seven whitepapers Bazavov et al. (2019); Brower et al. (2019); Cirigliano et al. (2019); Detmold et al. (2019); Joó et al. (2019); Lehner et al. (2019). This whitepaper covers the role of lattice QCD in neutrino-nucleus scattering, motivated principally by neutrino oscillations.

Neutrino-nucleus scattering experiments provide an abundance of information on neutrino masses and flavor mixing, on nucleon and nuclear structure, and on non-standard interactions between neutrinos and ordinary matter. To interpret these experiments cleanly, the key problem is to reconstruct the incident neutrino energy. The nuclear remnant is not, in these experiments, detected. It is therefore impossible to reconstruct the neutrino energy without modeling the nucleus in some way. This problem is complex, because it spans a range of energies—from hundreds of keV to a few GeV—that probe all aspects of the target nucleus.

The presence of many energy scales implies that a variety of theoretical techniques must work in concert. A convenient, organizational framework is nuclear many-body theory, which takes nucleonic properties as inputs. In this whitepaper, we discuss how these nucleonic properties can be obtained directly from the QCD Lagrangian using numerical simulations of lattice gauge theory. Although lattice QCD cannot settle every question in neutrino-nucleus scattering, it is reasonable to demand that our understanding of these processes be consistent with QCD. In many cases, the most straightforward route to the needed QCD knowledge is lattice QCD.

In this whitepaper, we discuss several calculations that should, as they mature, be incorporated into nuclear theory and neutrino event generators. A very important and very feasible example is the axial form factor of the nucleon. Lattice QCD has a notable history of calculating this and related observables, and calculations with full control of the systematic uncertainties are now coming of age. Here, “full control of systematic uncertainties” implies that a complete error budget is provided. The axial form factor is relatively straightforward: completely analogous calculations of vector form factors are possible with the same (indeed, overlapping) computational effort. The vector form factors have been measured in electron-proton and -neutron scattering, so an apt crosscheck is close at hand. Experience from form factors in meson physics suggests a simple, model-independent way to transmit the output of lattice QCD to event generators and, thus, analysis of experimental data.

Form factors of nucleons are only the beginning. Future oscillation experiments span beam energies such that computationally more demanding information is required. Just at the nucleon level, transition form factors to multibody final states are needed. For an inclusive data set, the object of interest is the nuclear hadron tensor, which can be obtained by combining the nucleonic hadron tensor from lattice QCD with a nuclear spectral function. In the deep inelastic region, new ways of computing parton distribution functions in lattice QCD are an exciting development. A further emerging component of lattice QCD consists of calculations of the properties of small nuclei—up to 4He today and to 6Li with exascale computing—can be used to test nuclear many-body theory and provide information via chiral effective theories to pin down the nuclear physics.

Lattice-QCD calculations with nucleon and nuclei are more challenging than the corresponding ones for mesons, because of unavoidable technical challenges that increase with the number of quark lines. Consequently, to perform the requisite calculations, improvements in methodology, algorithms, and software will be essential. Even assuming continuing ingenuity on those fronts, much of the work will require exascale computing resources. As in the past, a combination of high-capability and high-capacity computing will be needed. The former is needed for timely solution of mature problems, while the latter is necessary for developing feasible techniques for the challenging calculations, before making the jump to supercomputer centers.

I Introduction

Along with the first observation of the Higgs boson and the mounting evidence for dark matter, the discovery that neutrinos change flavor is one of the major advances in particle physics over the past twenty-five years. The discovery hinged on studies of neutrinos produced at the upper edge of the earth’s atmosphere Fukuda et al. (1998) and also explained a deficit in electron neutrinos from the sun Ahmad et al. (2001). These findings prompted an accelerator-based experimental program in Europe, Japan, and the United States, to make more accurate measurements of, for example, the squared mass differences. The increase in precision and sensitivity expected in future experiments raises the question whether the theoretical description of the relevant experiments must be further refined to exploit the new measurements to the fullest. In particular, as future, ambitious, long-baseline neutrino-oscillation experiments such as LBNF/DUNE Acciarri et al. (2015) and HyperK Abe et al. (2018) have come into focus, the quantification of uncertainties from the hadronic and nuclear physics of the detectors have become increasingly relevant. To this end, the lattice-QCD community has identified a set of feasible calculations that will be of special relevance. This program is described in this whitepaper.

An important goal of the experimental neutrino-physics program is to test the three-neutrino paradigm of the Standard Model. In this context, the Standard Model must be extended to allow for lepton flavor change. The simplest choice consistent with the standard gauge symmetries is to introduce a set of right-handed neutrino fields. Then lepton-flavor mixing and neutrino masses arise in the same way as in the quark sector, namely through Yukawa couplings to the Higgs field with a nonvanishing vacuum expectation value. To couple to the Higgs and left-handed-lepton doublets, the right-handed neutrino fields have to be gauge singlets. But then no symmetry principle forbids a mass term connecting neutrinos to themselves (i.e., of the kind first noted by Majorana Majorana (1937)), in contrast to the Higgs-generated Dirac mass term, which connects neutrino to antineutrino. The lack of direct evidence for right-handed neutrinos suggests that in this scenario the Majorana mass $M$ might be very large. If one supposes that the neutrino Yukawa couplings are not much different from light quarks or charged leptons, the propagating neutrinos have mass close to $M$ and to $m_{\nu}\approx y^{2}v^{2}/2M$, where $y$ is a Yukawa coupling and $v$ is the vacuum expectation value of the Higgs field. This mass hierarchy, known as the see-saw mechanism, provides a possible explanation of the tiny size of neutrino masses Minkowski (1977); *Yanagida:1980xy; *GellMann:1980vs. For example, if $M$ is a grand-unified mass scale around $10^{15}$ GeV, then $m_{\nu}\lesssim 0.03$ eV (for $y\lesssim 1$).

This theoretical framework means that the three-neutrino paradigm can be tested by measuring the neutrino mass-squared differences and the mixing angles and $CP$ violating phases of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix Pontecorvo (1957); *Pontecorvo:1967fh; *Maki:1962mu. Like the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix Cabibbo (1963); *Kobayashi:1973fv, the PMNS has three mixing angles. If the Majorana mass term appears, the PMNS matrix has three $CP$-violating phases instead of one as in the CKM matrix. The mixing angles and the CKM-like $CP$-violating phase can be measured in oscillation experiments, while the extra phases and the Majorana nature of neutrinos can be probed via the neutrinoless double-beta ($0\nu\beta\beta$) decay of certain nuclei. For lattice-QCD calculations relevant to $0\nu\beta\beta$, see the companion whitepaper “The Role of Lattice QCD in Searches for Violations of Fundamental Symmetries and Signals for New Physics” Cirigliano et al. (2019); here, the focus is on lattice-QCD research that will impact the oscillation experiments.

Oscillation experiments measure the energy spectrum of a neutrino beam after it has travelled a certain baseline distance. Unfortunately, neutrino beams have a wide energy spectrum, as shown in Fig. 1, so the center-of-mass energy of a collision is not known.

In contrast, quark-flavor experiments, for which lattice QCD has been crucial, study decays of strange, charmed, or $b$-flavored hadrons of precisely known mass. Here, the energy of the incident neutrino must be inferred from measurements of the final state. The targets in neutrino experiments are medium- to large-sized nuclei, such as 12C, 16O, or 40Ar, the remnants of which are not, in practice, be detected. That means that the mapping between final-state measurements and the initial energy inevitably requires theoretical knowledge of the neutrino interaction with the struck nucleus.

Consistency with QCD is a clearly desirable characteristic of nuclear models used to deduce the connection between final and initial states. Thus, it makes sense to incorporate lattice QCD as soon as results with full, reliable error budgets are available. As discussed in more detail in Ref. Alvarez-Ruso et al. (2018), the nuclear models rely in part on properties of the nucleon as inputs. Many of these quantities can be calculated in lattice QCD in the near term, with the precision depending on the quantity. Of course, single-nucleon calculations are not in themselves enough. Calculations of the properties of multi-nucleon systems must be developed concurrently and, once mature, also incorporated into the nuclear modeling.

The theory behind neutrino-nucleus collisions is complex because it spans a range of energies that probe all aspects of the target nucleus. Nuclear excitation energies are, typically, dozens of keV, while the average binding energy is 8.6 MeV (in ${}^{40}\text{Ar}$), and the typical Fermi motion of a nucleon is around 250 MeV. In the regime relevant to oscillation experiments, the energy transfer to the nucleus ranges between $\sim$200 MeV and the neutrino energy itself, although much of transferred energy is carried off by nucleons and pions, rather than the nuclear remnant. Thus, it is a challenge to arrive at a comprehensive approach to the entire problem. Most approaches start with nuclear many-body theories, in which the nucleus is described by a nuclear wave function of a collection of interacting nucleons; see, for example, Ref. Benhar and Lovato (2015); Carlson et al. (2015). It is at this point in the analysis that nucleon-level matrix elements enter. One should bear in mind, however, that single-nucleon physics is not enough: multi-body effects are needed for scattering events that knock out two (or more) nucleons. Even in nuclear spectroscopy, three-body potentials improve the agreement with observed nuclear levels Epelbaum et al. (2009); Machleidt and Entem (2011); Carlson et al. (2015). Often these calculations use phenomenological potentials, but effective field theory (EFT) offers a direct connection to QCD Weinberg (1990); van Kolck (1994); Kaplan et al. (1998); Meißner (2016). Chiral EFTs are, however, limited to a kinematic range where the momenta are small relative to the chiral symmetry breaking scale $\Lambda_{\chi}\sim 700~{}\text{MeV}$. Even then, the reliability of the application of nuclear EFT to large atomic number systems, such as argon, requires significant development, testing, and, eventually, verification.

Juszczak et al. (2006); *Zmuda:2015twa

Andreopoulos et al. (2010); *Alam:2015nkk These issues are further intertwined with the constraints of how event generators Juszczak et al. (2006); Hayato (2009); Casper (2002); Andreopoulos et al. (2010); Gallmeister et al. (2016) and detector simulations are implemented. Inconsistencies arise in the current approach where, for example, the axial form factor of the nucleon is extracted from $\nu A$ scattering data assuming one nuclear model and then used in event generators employing another.

A central goal of nuclear theory in this arena should therefore be to define a path forward that allows for a quantified nuclear uncertainty to be presented for experiments such as DUNE and HyperK. Achieving this is a challenging task and will require input and constraints from lattice QCD in order for it to be successful. In addition to the single- and few-nucleon amplitudes noted above, it will be valuable to compute directly the properties of small nuclei. At present, calculations involving nuclei up to 4He are possible. In addition to being interesting in their own right, such lattice-QCD calculations of few nucleon systems can be used to constrain low energy constants (LECs) in the EFTs. This approach has already been applied to static quantities, such as magnetic moments. A next step will be to work with matrix elements of electroweak currents, to build up effects associated with two- and higher-body contributions, as well as more complex contributions such as pion production. In combination with experimental constraints from $eA$ scattering, and neutrino scattering on light nuclear targets,111Indeed, recent discussions of future experiments with deuterium or hydrogen targets Morfín (2018) hinge on noting the utility of nucleon-level amplitudes in nuclear many-body theory. it is hoped a robust uncertainty can be determined.

To study neutrino oscillations, we are interested in the processes

[TABLE]

where $A$ denotes the nucleus and $X$ the combination of all final-state hadrons including the remnant of the nucleus. The charged weak current responsible for these interactions has the well-known $V-A$ structure. Properties of the vector current can be inferred from electromagnetic scattering, up to isospin corrections (which are negligible for the needed precision; see Sec. IV). On the other hand, because the weak charge of the proton is so small, $Q_{\text{w}}^{p}=0.0719\pm 0.0045$ Androić et al. (2018), at the energies of interest, only neutron-neutrino (and proton-antineutrino) scattering is sensitive to the axial current. These circumstances offer the possibility of testing lattice-QCD methodology with the vector current before relying on it for the axial current.

The quantity needed to describe the strong-interaction side of the scattering depends on the energy transferred. At the lowest energies, the only possibility is coherent elastic scattering via the weak neutral current, with $X=A$ Freedman (1974); Brice et al. (2014). Coherent neutrino-nucleus interactions have recently been observed for the first time Akimov et al. (2017). As the energy increases slightly, the excitation spectrum of $A$ is traced out: $X=A^{*}$. The needed quantities are matrix elements between different nuclear levels. In lattice QCD, one would have to simulate the whole nucleus directly, which is currently feasible only for nuclei much smaller than those in the cesium-iodide detector of Ref. Akimov et al. (2017).

At high enough (but still low) energy, a single nucleon can be knocked out. At its heart, the scattering is

[TABLE]

with the initial and final-state nucleons in the nuclear environment. Such scattering off of a constituent in a bound-state without extra particles is known as quasielastic. Then nuclear many-body theory requires single-nucleon matrix elements of the form $\langle p(p^{\prime})|J_{\nu}|n(p)\rangle$, between a neutron of momentum $p$ and a proton of momentum $p^{\prime}$ (or the $p\to n$ counterpart for antineutrino beams). These matrix elements are straightforward to calculate in lattice QCD; see Sec. II. If pions can be produced, the final state can be a $\Delta(1232)$ resonance, an excited nucleon $N^{*}$, or a two-body state $N\pi$. In the experiment, these all end up as $N\pi$ so their amplitudes interfere. In fact, lattice QCD can provide not only the associated transition matrix elements, in the idealization of the resonance as a stable particle (e.g., $\langle\Delta^{+}|J_{\nu}|n\rangle$), but also enough information to describe the full multi-hadron nature of the final state (at least up to further inelasticities); see Sec. III. The quasielastic and resonance regions overlap, because the kinetic energy of Fermi motion is a bit larger than the pion mass. This overlap is illustrated with experimental data in Fig. 2.

Another contribution in this region arises from many-body nuclear dynamics, for example, when the probe interacts with pairs of correlated nucleons. This contribution is described by “two-body currents” (see Refs. Carlson and Schiavilla (1998); Bacca and Pastore (2014) and references therein). Now a further set of matrix elements is needed, namely of the form $\langle NN|J_{\nu}|NN\rangle$. Note that in QCD language, the same current is employed just for two-body initial and final states.

Once the energy is high enough to produce several pions, it is not possible to enumerate every final-state hadron. In this case, however, lattice QCD can be used to compute nucleon and nuclear structure functions. In high-energy physics, structure functions are most familiar in deep-inelastic scattering, where the operator-product expansion (OPE) can be used. Lattice-QCD calculations can be used to determine the moments of the parton distribution functions (PDFs) that enter in the deep-inelastic region, and indeed extraction of the full dependence of PDFs on the longitudinal momentum fraction, $x$, is becoming possible Lin et al. (2018). Moreover, the definition of structure functions is very general. Lattice QCD may be an ideal way to compute them in the so-called shallow-inelastic region with energy above the resonance region but insufficient for the OPE; see Sec. III.

In summary, then, the goals for lattice QCD for neutrino oscillation physics are to calculate matrix elements of the form

[TABLE]

where the initial and final states are single nucleons, two nucleons, nucleons with a pion (including resonances), or small nuclei. In the last case, $\mathcal{O}$ denotes an operator appearing in the OPE, or a bilocal, spatially-separated operator arising in the calculation of PDFs. The lattice-QCD calculations of these and related matrix elements have a long history, motivated principally by the desire to understand nucleon and nuclear structure. For a broad survey, see our companion whitepaper “Hadrons and Nuclei” Detmold et al. (2019).

Bazavov et al. (2013); *Bazavov:2017lyh Recall that lattice QCD calculates hadronic correlation functions, which contain information about the masses and matrix elements of interest; the information is extracted by fitting the behavior of the correlation functions in (Euclidean) time. Several technical difficulties make baryon calculations more difficult than the corresponding calculations for mesons. First, statistical errors on baryon correlation functions are larger and more poorly behaved in time Parisi (1984); Lepage (1989); Wagman (2017). Second, it has proven more difficult, in practice, to disentangle matrix elements of the ground-state baryons from that of their excitations Owen et al. (2013). Last, the dependence of baryon properties on the light quark mass (used in the simulation) is less well described by the low-energy EFT of pions and baryons. All these difficulties can be addressed with more computing. The signal-to-noise problem can clearly be attacked with higher statistics. It can also be mitigated by choosing more sophisticated operators to create and annihilate baryon states; this method is also the way to better filter out the excited states. Finally, more computing also enables simulations with lighter and even physical quark masses Dürr et al. (2011); Bazavov et al. (2013); Blum et al. (2016).

The rest of this whitepaper is organized as follows. In Sec. II, we discuss calculations that are relatively straightforward. These include nucleon form factors, which are needed to describe quasielastic scattering, and moments of PDFs, which are needed in the deep-inelastic region. We discuss the form factors in considerable detail, because the time to incorporate these results into event generators is soon or, arguably, now. In particular, having the correct slopes for the form factors is crucial to gaining quantitative control of the cross section. More challenging calculations are covered in Sec. III. This class of problems is large and varied: transitions to resonances and multibody states, calculations for shallow- and deep-inelastic scattering, and the vector and axial matrix elements of small nuclei. Section IV turns to calculations that are far enough beyond that state of the art that new ideas or computing facilities greater than exascale are needed. Foreseeable computing needs are covered in Sec. V, noting the separate needs for both capability and capacity computing.

II Straightforward calculations

The most straightforward matrix elements to calculate are those with one stable hadron in the initial state, and one or none in the final state. Here we focus on the matrix elements of electroweak currents, $\langle N|J_{\mu}|N\rangle$, which directly enter neutrino-nucleon scattering, and matrix elements of local operators, $\langle N|\mathcal{O}|N\rangle$, where $\mathcal{O}$ appears in the operator-product expansion of two $J$ currents, which arise in the analysis of deep-inelastic scattering.

II.1 Nucleon form factors

As discussed in Sec. I, neutrino-nucleon scattering, Eq. (2) is a key process even though the target is a nucleus. The $V-A$ charged current of interest is $J_{\mu}^{+}=\bar{u}\gamma_{\mu}(1-\gamma_{5})d$. The matrix element for $n\to p$ can be decomposed into Lorentz covariant combinations of momentum and spin, multiplied by form factors Bhattacharya et al. (2011):

[TABLE]

where $M_{N}=(M_{p}+M_{n})/2$, $q=p^{\prime}-p$ and $\bar{u}$ and $u$ are associated spinor factors. $F_{1}^{\text{CC}}(q^{2})$, $F_{2}^{\text{CC}}(q^{2})$, $F_{A}^{\text{CC}}(q^{2})$, and $F_{P}^{\text{CC}}(q^{2})$ are known as the Dirac, Pauli, axial, and induced pseudoscalar form factors, respectively. The induced scalar and tensor form factors, $F_{S}^{\text{CC}}(q^{2})$ and $F_{T}^{\text{CC}}(q^{2})$, are suppressed by $G$ parity violation; they are known as second-class currents Weinberg (1958). For neutral-current processes, additional form factors $F_{i}^{\text{EM},N}$ and $F_{i}^{\text{NC},N}$ are needed: the charged-currents are all isovector, but the neutral currents contain an isoscalar contribution as well. Here, $N$ denotes either a proton $p$ or neutron $n$.

Because the up- and down-quark masses are so similar, isospin violation can be neglected and, thus, the charged-current form factors of the vector current (i.e., Dirac and Pauli) can be related to their electromagnetic counterparts, up to small corrections from isospin violation. The Dirac and Pauli form factors are usually re-expressed as electric, $G_{E}(q^{2})=F_{1}(q^{2})+q^{2}F_{2}(q^{2})/(M_{n}+M_{p})^{2}$, and magnetic, $G_{M}(q^{2})=F_{1}(q^{2})+F_{2}(q^{2})$, form factors (even for CC and NC). Expressions relating the differential neutrino-nucleon cross section to the form factors can be found, for example, in Refs. Llewellyn Smith (1972); Formaggio and Zeller (2012).

Most neutrino scattering experiments are performed in a kinematic region of a few GeV, so tracing out the full $q^{2}$ dependence is possible and desirable (see below). Below 1 GeV it is convenient to focus attention on the intercepts $F_{i}(0)$ and (conventionally normalized) slopes

[TABLE]

for $i\in\{A,S,T{,P}\}$. The quantities $r_{i}$ are usually called “radii”, although the neutron’s $r_{E}^{2}$ is negative.

A precise knowledge of the charged-current versions of these quantities is essential for determining the neutrino-nucleon cross section. The intercepts and slopes of $G_{E}^{\text{CC}}$ and $G_{M}^{\text{CC}}$ are well determined from electromagnetic processes and isospin relations. Further, the intercept $F_{A}^{\text{CC}}(0)=g_{A}=-1.2723(23)$ is known from neutron $\beta$ decay Tanabashi et al. (2018). The axial coupling $g_{A}$ has been calculated in lattice QCD, although it will be some time before it can be computed with comparable precision to experiment. Nevertheless, it is an extremely important benchmark, and once the lattice-QCD precision becomes competitive with experiment, the result could clear up some puzzles surrounding neutron-decay measurements (see below).

On the other hand, the axial-charge radius-squared $r_{A}^{2}$ is less well known. Historically, the axial form factor has been fit to the so-called “dipole” form:

[TABLE]

such that $r_{A}^{2}=12/m_{A}^{2}$. Experiments report this “axial mass”, $m_{A}$, so a comparison of reported values illustrates the current status. It has been extracted from quasielastic scattering on deuterium targets, finding (e.g.) $m_{A}=1.02(3)$ GeV Bodek et al. (2008), and from pion electroproduction, finding $m_{A}=1.08(4)$ GeV Liesenfeld et al. (1999); Bernard et al. (1992). More recent experiments find larger values: $m_{A}=1.20(12)$ GeV at K2K Gran et al. (2006), $m_{A}=1.27(15)$ GeV at MINOS Dorman (2009), and even $m_{A}=1.35(17)$ GeV at MiniBooNE Aguilar-Arevalo et al. (2010a), in neutrino charged-current quasielastic scattering with water, iron, and mineral-oil targets, respectively. With 2p-2h corrections, however, NOMAD Lyubushkin et al. (2009), with a Kevlar target, finds $m_{A}=1.05(6)$ GeV and

MINERvA Fields et al. (2013); *Fiorentini:2013ezn, with a carbon target, finds the quasielastic cross section to be compatible with $m_{A}=0.99$ GeV. Note that all of these determinations of $m_{A}$ assume a nuclear model for the target material, which is not the same among the various collaborations. Moreover, nuclear modeling uncertainties typically come only from varying parameters of their choice model, not from studying comparisons among different models.

The uneasy agreement of these results can be removed by switching to a model-independent parametrization of $F_{A}(q^{2})$ Bhattacharya et al. (2011). For example, a reanalysis of 1980s deuterium bubble-chamber data Meyer et al. (2016a) finds $\sqrt{12/r_{A}^{2}}=1.01(24)~{}\text{GeV}$. These data are chosen because the nuclear model of the deuteron is under relatively good control. The main conclusion of Ref. Meyer et al. (2016a) is that introducing only one free parameter with a qualitatively acceptable but conceptually incorrect shape, as in Eq. (6), leads to gross underestimates of the uncertainty, even when the fit quality is high.

Figure 4, from Ref. Hill et al. (2018), shows the dependency of $\nu n$ quasielastic cross section on $E_{\nu}$, assuming $r_{A}^{2}$ is known with $20$% uncertainty. As one can see, this quantity affects the both the normalization and fall-off of the cross section, which are needed, respectively, to determine the mixing angle and mass difference in an oscillation. Furthermore, a lattice-QCD calculation with 20% uncertainty (compared to 50% in Ref. Meyer et al. (2016a)) is an important milestone, because then the $r_{A}^{2}$ uncertainty becomes subdominant, at least until other uncertainties have been reduced.

The lattice-QCD community has been pursuing the calculation of the nucleon form factors for a long time. A representative set of recent work can be found in Refs. Bali et al. (2015); Bhattacharya et al. (2016); Green et al. (2017); Alexandrou et al. (2017a); Capitani et al. (2019); Gupta et al. (2017); Jang et al. (2018a); Yamanaka et al. (2018); Chang et al. (2018a); Ishikawa et al. (2018); Liang et al. (2018a); Ottnad et al. (2018); Bali et al. (2019); Shintani et al. (2019); Constantinou et al. (2018) Significant improvements have been made to investigate the quark-mass, finite-volume, and finite-lattice-spacing dependence, and the effects of excited-state contamination in the correlation functions. With these technical and algorithmic advances, lattice QCD can calculate not only the isovector contribution but also the computationally more demanding isoscalar and strange-quark contributions, which are needed for neutral-current processes, discussed below.

Sample lattice-QCD calculations Jang et al. (2018a, b) of the nucleon isovector electric and axial form factors—$G_{E}$ and $F_{A}$—are shown in Fig. 4.

Eight different $2+1+1$-flavor HISQ ensembles generated by the MILC collaboration Bazavov et al. (2013); *Bazavov:2017lyh with lattice spacings in the range 0.06–0.12 fm and pion mass in the range 130–310 MeV are employed. In this calculation, excited-state contamination is controlled via a three-state fit. The results are in good agreement with the experimental data for the nucleon electromagnetic form factor $G_{E}(q^{2})$ On the other hand, the axial form factor is not as steep as experimental determinations with $m_{A}\approx 1$ GeV Bernard et al. (2002), yet is compatible with MiniBooNE’s $m_{A}\approx 1.35$ GeV Aguilar-Arevalo et al. (2010a). Despite the many laudable aspects of Ref. Jang et al. (2018a), a full and robust accounting of all systematics involved in these lattice-QCD calculations has not yet been feasible. Reliable confrontation with precise experimental data for $G_{E}$—and, hence, a solid prediction of $F_{A}$—requires an increase in computational resources to overcome the technical obstacles to nucleon matrix elements, discussed in Sec. I.

The status of lattice-QCD calculations of $g_{A}$ and $r_{A}^{2}$ is shown in Fig. 5.

The left plot Gupta et al. (2018), for $g_{A}$, shows that lattice-QCD is at this time much less precise than the results from neutron $\beta$ decay.222The color code here is adapted from the Flavor Lattice Averaging Group Aoki et al. (2017); *Aoki:2019cca; *FLAGweb, as specified in the Appendix of Ref. Bhattacharya et al. (2015). Note, however, that bottle and beam experiments measuring the neutron lifetime yield values of $g_{A}$ that differ by $3\sigma$. For example, a 2015 bottle measurement leads to $g_{A}=1.2749(11)$ Arzumanov et al. (2015), while a 2013 beam measurement leads to $g_{A}=1.2684(20)$ Yue et al. (2013). It would be interesting to know the answer from lattice QCD. The precision required depends on whether the (average of several) calculation(s) lands between the two neutron-lifetime values or outside the interval. In the latter case, at least percent-level precision is needed, which is likely to be achieved with three years (assuming sustained computing support). If lattice QCD lands in the middle, 0.3% precision is needed. In this scenario, we would also need 1+1+1(+1)-flavor ensembles, since the isospin symmetry would play an important role at such precision; it would take 5–10 years to account for full systematics.333Note that the normalization of the matrix element can be blinded with an multiplicative offset Meyer et al. (2016b), to guard against analyst bias. The results in Fig. 5 (left) have not, however, employed this technique.

The right plot Hill et al. (2018), for $r_{A}^{2}$, shows significant problems: the analysis with the $z$ expansion Meyer et al. (2016a) debunks the uncertainty estimates of determinations predicated on the dipole form. The model independent results (red; between the horizontal lines) illustrate the best estimate of $r_{A}^{2}$ without such strong assumptions. One should bear in mind that the “experimental” determinations all make assumptions: without new $\nu d$ and $\bar{\nu}p$ experiments Morfín (2018), it seems nearly impossible to improve the situation via experiment. On the other hand, lattice gauge theory can provide an ab initio result from QCD. Indeed, lattice QCD is beginning to play a role, but another generation of calculations is needed before fully definitive results with uncertainties small enough to make an impact on cross section calculations are achieved.

For the full energy range of LBNF/DUNE, it will be necessary to trace out the full $q^{2}$ dependence of the form factors. It is imperative to use a model-independent parametrization based on general analytic properties. In the complex-$q^{2}$ plane, the vector (axial) form factors have a cut starting at $q^{2}=t_{\text{cut}}\equiv 4M_{\pi}^{2}$ ($q^{2}=t_{\text{cut}}\equiv 9M_{\pi}^{2}$) and extending to $\infty$ on the real axis. The cut lies outside scattering kinematics $q^{2}<0$ but nevertheless prevents a useful series expansion in $q^{2}$ around the origin. A rigorous way to proceed is to introduce a conformal mapping that maps the cut to the unit circle Meiman (1963); Hill and Paz (2010):

[TABLE]

where the parameter $t_{0}$ can be chosen to center the $q^{2}$ range of interest on $z=0$; in general, spacelike $q^{2}\to-\infty$ maps to $z\to 1$. An expansion of the form

[TABLE]

thus has an expansion parameter $|z|<1$. Moreover, unitarity in quantum mechanics ensures that the series is uniformly convergent on this interval. In fact, unitarity leads to bounds on the coefficients $a_{k}$ that the dipole form, Eq. (6), violates Bhattacharya et al. (2011).

In practice Bhattacharya et al. (2011); Hill (2006); Bernard et al. (2009); Bailey et al. (2009), the $z$ expansion converges after a few terms. Because on the nonlinear mapping, even an intercept and slope in $z$ give a form factor with a physical shape (i.e., similar to those shown in Fig. 4). As lattice data improve, more and more terms will become resolvable. As in CKM physics Bernard et al. (2009); Bailey et al. (2009), lattice-QCD papers can provide the coefficients, their uncertainties, and their correlations; several lattice-QCD calculations of $F_{A}$ do the same Alexandrou et al. (2017a); Capitani et al. (2019); Green et al. (2017); Gupta et al. (2017). Finally, code for taking such $z$-expansion input is included in the GENIE event generator Andreopoulos et al. (2010) module for the axial form factor, and work is underway to extend this to the vector form factor channel.

Although not crucial to neutrino oscillations, the same experiments study weak neutral-current interactions of the $Z$ boson and, possibly, non-Standard bosons de Gouvêa and Kelly (2016); Heeck et al. (2018). The corresponding Dirac and Pauli form factors can be obtained from the proton and neutron electromagnetic form factors and the strange-quark contribution (accessible in parity-violating elastic electron-scattering experiments Armstrong and McKeown (2012)) as

[TABLE]

Using the most recent $z$-expansion fit to nucleon electromagnetic form factors Ye et al. (2018) and a new lattice-QCD calculation of strange-quark form factors Sufian (2017), one can see that the strange-quark contribution increases the neutral-current Pauli form factor, $F_{2}^{\text{NC}}(q^{2})$, by about 3.1% and 2.5% at $q^{2}=0$ and $q^{2}=-0.1~{}\text{GeV}^{2}$, respectively. Although the strange-quark contribution is small, the coefficients $(\frac{1}{2}-\sin^{2}\theta_{\text{W}})$ and $\sin^{2}\theta_{\text{W}}$ suppress the two combinations of nucleon electromagnetic form factors in Eq. (9), such that the strange-quark sea makes an important contribution to $F_{2}^{\text{NC}}(q^{2})$ at low $q^{2}$.

Similarly, assuming isospin symmetry and the absence of second-class currents, one can relate the neutral-current axial form factor to the charged-current axial and strange-quark axial form factors Garvey et al. (1993a, b):

[TABLE]

It has been shown Aguilar-Arevalo et al. (2010b, 2015) that the effect of Pauli blocking becomes very significant in the region $0<-q^{2}\lesssim 0.2~{}\text{GeV}^{2}$. Therefore, a precise lattice-QCD calculation of $F^{\text{NC}}_{A}(q^{2})$ is required for a precise estimate of the neutral-current (anti)neutrino-nucleon scattering cross section.

Finally, we note that quasielastic neutrino and antineutrino scattering would be sensitive to the presence of the second-class currents, $F_{S}$ and $F_{T}$ in Eq. (4), characterized by a different $G$-parity to the standard vector and axial currents of the Standard Model. The search for such currents has long been pursued in the $\beta$-decay experiments and in muon-capture experiments, but the measurement of polarization observables in the quasielastic scattering both of nucleons and of hyperons has been shown to be sensitive both to $G$ invariance and to $T$-invariance Fatima et al. (2018). Lattice QCD can contribute to these tests through calculations of induced scalar and tensor currents, including calculations of transition form factors to the rest of the SU(3) baryon octet ($\Lambda$ and $\Sigma$ as well as $p$ and $n$), such as those in Refs. Lin (2009); Sasaki (2014).

II.2 Moments of parton density functions

Lattice QCD can be used to calculate matrix elements of other operators besides the electroweak currents. An important class of operators are those the appear in the operator-product expansion of two currents. Their matrix elements are related to the moments of structure functions in deep-inelastic scattering. For a full discussion, see the USQCD companion white paper “Hadrons and Nuclei” Detmold et al. (2019). Here, applications to neutrino physics are discussed.

In 2001, the NuTeV collaboration determined the on-shell weak mixing angle, $\sin^{2}\theta_{W}\equiv 1-m_{W}^{2}/m_{Z}^{2}$, to be $0.2277\pm 0.0013_{\text{stat}}\pm 0.0009_{\text{syst}}$ Zeller et al. (2002) in deep-inelastic neutrino scattering off iron. This result is $2.7\sigma$ discrepant from the current world average of other experiments, $0.22343\pm 0.00007$ Tanabashi et al. (2018). This discrepancy, which is known as the “NuTeV anomaly”, has no universally accepted explanation, although many possibilities have been raised Londergan and Thomas (2003); Ding et al. (2005); Gluck et al. (2005); Eskola and Paukkunen (2006); Cloet et al. (2009); Bentz et al. (2010).

One suggestion that may account for part of the anomaly is the strange-antistrange parton asymmetry Davidson et al. (2002); Kretzer et al. (2004), $\langle x\rangle_{s_{-}}=\int dx\,x\left[s(x)-\bar{s}(x)\right]$, where $s(x)$ ($\bar{s}(x)$) is the (anti)strange parton distribution function, as a function of parton momentum fraction $x$. A global analysis of several experimental data sets gives $\langle x\rangle_{s_{-}}\approx 0.0018$ Lai et al. (2007), which is consistent with a 2006 NuTeV analysis of dimuon production Mason (2006). The global analysis does not, however, find a tight constraint: the authors of Ref. Lai et al. (2007) present the range $-0.001<\langle x\rangle_{s_{-}}<0.005$ at 90% confidence level.

In view of the uncertain of $\langle x\rangle_{s_{-}}$ from global fitting, a first-principles lattice-QCD calculation is warranted. There is, however, no local operator which corresponds to $\langle x\rangle_{s_{-}}$. Instead, one can calculate the third moment from the local operator $\bar{s}\gamma_{\mu}D_{\nu}D_{\lambda}s$ which corresponds to $\langle x^{2}\rangle_{s_{-}}=\int dx\,x^{2}(s(x)-\bar{s}(x))$. Assuming $s(x)-\bar{s}(x)$ changes sign only once, $\langle x^{2}\rangle_{s_{-}}$ should give the same sign as that of $\langle x\rangle_{s_{-}}$. This quantity can also be used to constrain the $x$-dependent distribution, but since it is expected to be small, calculations will require significant resources.

III Challenging calculations

In this section, calculations that are computationally more difficult than the form factors in Sec. II are discussed. That said, the conceptual formalism underlying these calculations is well established, and pilot calculations provide some idea of how more complete calculations can be carried out. More complicated final states in the resonance regions (Sec. III.1), the shallow inelastic region (Sec. III.2), and the deep inelastic region (Sec. III.3) are discussed, as are calculations of the axial charge, and related quantities, of small nuclei (Sec. III.4).

III.1 Transition form factors: resonances and multibody final states

Neutrino scattering above the pion-production threshold constitutes the resonance region, where the scattered nucleon is excited into resonances, beginning with the $\Delta(1232)$. To describe the data in this regime thus requires a quantitative knowledge of the $N\to\Delta$ and $N\to N^{*}$ transitions, mediated through an external current. Because these hadrons are unstable, they can also be viewed as a nucleon with one or more pions, which are the only hadrons composed of the light $u/d$ quarks stable under the strong interaction.

Lattice QCD has a long history of calculations of the transition form factors to the $\Delta$, treating it as stable. Both the vector current Leinweber et al. (1993); Alexandrou et al. (2008), and the axial current Alexandrou et al. (2007) have been studied with unphysically large quark masses, such that $M_{\Delta}$ at these quark masses lies below the $N\pi$ threshold. These calculations are useful benchmarks for comparisons with non-lattice approaches that neglect the two-body nature of the resonance. Although not as rigorous as the methods discussed below, this “quick and dirty” approach may be timely, for example, providing qualitative input to understand better the MiniBooNE low-energy backgrounds from $\Delta\to N\gamma$ Aguilar-Arevalo et al. (2018).

Because of the finite volume and Euclidean signature, calculations with two-body states in lattice QCD are conceptually and computationally more difficult Lüscher (1986); *Luscher:1990ux; Lellouch and Lüscher (2001) than the calculations discussed in Sec. II. For example, the Lüscher method Lüscher (1986); *Luscher:1990ux relating energy shifts at finite volume to infinite-volume momentum-dependent phase shifts has been used to study the $\rho$ meson Aoki et al. (2007); Feng et al. (2011); Dudek et al. (2013); Wilson et al. (2015); Alexandrou et al. (2017b), as well as $I=2$ $\pi\pi$ phase shifts Beane et al. (2008a); Feng et al. (2010); Beane et al. (2012); Dudek et al. (2011, 2012) from first principles. The theoretical framework for understanding the transition to multihadron states from Euclidean-space lattice QCD calculations have been further developed over the past several years. Notably, the formalism has been extended both to inelastic scattering Detmold and Savage (2004); He et al. (2005); Hansen and Sharpe (2012); Briceño and Davoudi (2013a); Guo et al. (2013) with several two-body channels, and to three-body scattering Kreuzer and Hammer (2011); Briceño and Davoudi (2013b); Meißner et al. (2015); Briceño et al. (2017, 2019); Döring et al. (2018), and there have now been several computational applications of these advances Beane et al. (2008b); Dudek et al. (2016); Briceño et al. (2018a); Woss et al. (2018).

A quantitative understanding of resonance production entails extending the formalism to encompassing transitions mediated through external currents, corresponding here to both vector and axial currents. The needed formalism to two-body final states, and for arbitrary spin, has now been developed Briceño and Hansen (2015). The applications have largely focused on the meson sector. To cite an example bearing some similarity to $W^{*}n\to\Delta$ in neutrino scattering, the $\gamma^{*}\pi\to\rho$ transition has been computed in lattice QCD Briceño et al. (2016); Alexandrou et al. (2018), providing the first rigorous calculation of the transition form factor to an unstable hadron, illustrated in Fig. 6

In addition, methods to extract resonance-to-resonance transitions, for example, $\gamma^{*}\rho\to\rho$, via lattice calculations of two-to-two transition amplitudes, in this case $\gamma^{*}\pi\pi\to\pi\pi$, have been developed Briceño and Hansen (2016); Baroni et al. (2018). This opens the possibility for calculations of two-body currents, that is, matrix elements of the form $\langle NN|J^{\mu}|NN\rangle$ needed for two-nucleon knockout.

Thus, the theoretical underpinnings for understanding resonance production in $\nu N\to\Delta$, $N^{*}$, and $N\pi$, are therefore largely in place. Calculations of multihadron states containing baryons are complicated by the extra complexity of the systems relating to the increased number of quarks, by poorer signal-to-noise ratios, and by the larger number of open channels. Even so, the first ab initio determination of $\Delta(1232)$ resonance parameters appeared in 2017 Andersen et al. (2018), albeit for a simulation with quark masses corresponding to a pion mass of $280~{}\text{MeV}$, yielding a $\Delta$-$N$-$\pi$ coupling in agreement with phenomenological determinations. As the invariant mass of the system increases within the resonance regime and the pion mass is decreased to its physical value, however, inelastic processes and three- and higher-body final states become relevant. Further theoretical work is needed to encompass transitions to three or more particles, and further development of efficient algorithms is needed to evaluate the larger number of Wick contractions that increasingly dominate the computational cost of the calculation. The interplay between theoretical methods and practical algorithms is, of course, ideally researched on high-capacity computing facilities.

III.2 Hadron tensor for shallow and deep inelastic scattering

At higher energies, more and more pions are produced and a full theoretical description of any given final state becomes impractical. One can, however, study the sum over all final states via the optical theorem and consider forward matrix elements of the product of two currents. Whereas nuclear many-body theory decomposes the low-multiplicity cases into products of nuclear wavefunctions and nucleon (and $N\pi$, …) form factors, here one can decompose the nuclear hadron tensor, $\langle A|J_{\mu}^{\dagger}J_{\nu}|A\rangle$, into a spectral function Rocco et al. (2016) and the nucleon hadron tensor, $\langle N|J_{\mu}^{\dagger}J_{\nu}|N\rangle$. A recent development in lattice QCD is to calculate this quantity from a four-point correlation function. This approach is especially appealing in the region, sometimes called shallow inelastic, between the resonances and deep-inelastic scattering, where no other theoretical tool holds much promise Alvarez-Ruso et al. (2018).

The Euclidean hadronic tensor Liu and Dong (1994); Liu et al. (1999); Liu (2000); Aglietti et al. (1998); Detmold and Lin (2006); Liu (2016); Chambers et al. (2017); Hansen et al. (2017) can be decomposed in terms of structure functions that are related to their Minkowski counterparts through a Laplace transform. Thus, to obtain the desired structure functions, an inverse Laplace transform is needed, an ill-posed problem that arises in many fields. Three approaches to this problem are the maximum entropy method (MEM) Liu (2000), the MEM with a prior to stabilize the fit for Bayesian reconstruction (MEM-BR) Liu (2016), and the Backus-Gilbert method Hansen et al. (2017). These three numerical approaches have been studied recently Liang et al. (2018b). The Backus-Gilbert method yields a single broad peak in the energy spectrum from lattice data with 20 points in Euclidean time. With both the MEM and MEM-BR, the elastic peak and the resonance peak are resolved, with the MEM-BR producing sharper peaks and a more stable reconstruction. Given the test lattice spacing of 0.12 fm, there is no excitation spectrum above 2 GeV and, thus, no strength in the spectral weight above the resonance region. For a finer lattice with spacing 0.04 fm, the spectral weight up to 5 GeV is accessible. Even though it may still not be sufficient to resolve the individual resonances, the fact it can cover both the resonance and the shallow inelastic scattering regions makes the lattice hadronic tensor calculation a promising theoretical tool to address the total cross section of the neutrino-nucleus scattering over a wide range of energy transfers up to 5 GeV.

The hadron tensor can also be computed for deep-inelastic scattering. In this case, the calculation needs to be able to access the kinematic region where $Q^{2}>4~{}\text{GeV}^{2}$ and energy transfer $\nu>5$ GeV where the higher-twist contributions are suppressed. The Euclidean correlation function can also be analyzed with the OPE, along the lines of the suggestion for calculating the shape function of the inclusive $B$-meson semileptonic decay rate Aglietti et al. (1998). In addition, using a fictitious heavy-quark propagator between the currents to calculate moments has been proposed Detmold and Lin (2006). A related approach is also discussed in Ref. Chambers et al. (2017). Unlike the approaches discussed in the next subsection, the hadron-tensor approach to deep-inelastic scattering does not need to match to the infinite-momentum frame.

III.3 Parton densities for neutrino deep-inelastic scattering

The parton distribution functions (PDFs) will be important inputs in the upcoming precision neutrino-physics experiments, particularly at large Bjorken $x$ and at the highest energies of the DUNE beam, $\sim 4$–5 GeV. For these kinematics, current global-fit PDFs either suffer greatly from the theoretical uncertainty in their nuclear treatment or rely mainly on extrapolation from intermediate $x$.

Direct calculation of the Bjorken $x$ dependence of hadron structure in lattice QCD has only recently become possible thanks to the development of the large-momentum effective theory (LaMET) Ji (2013), which introduces a large momentum $P$ to connect Euclidean lattice QCD to the desired Minkowski distributions.444See Refs. Rossi and Testa (2018); Cichy and Constantinou (2018); Monahan (2018) for critical discussions of this approach. This framework has allowed the first direct lattice-QCD computations of the $x$ dependence of parton distributions Lin et al. (2015). Further developments spurred on by these developments include that of pseudo-PDFs Radyushkin (2017), and that of matrix elements of gauge-invariant current-current correlators Ma and Qiu (2018). In these new approaches, valence- and sea-quark structure can be disentangled, which leads to the possibility of using lattice-QCD calculations to directly compare with experiments on large-$x$ structure in SoLID at Jefferson Lab, with sea structure in Drell-Yan experiments at Fermilab or with data from a future electron-ion collider. In addition to the hadron-tensor method described in the previous subsection, these new approaches and numerical investigations thereof are described in detail in the companion whitepaper “Hadrons and Nuclei” Detmold et al. (2019).

The lattice-QCD effort so far has focused on isovector combinations of PDFs, that is, the difference between the up and down distributions. A recent joint lattice-QCD and global-fitting community report Lin et al. (2018), an effort led by USQCD members, demonstrated that a calculation of the isovector proton PDF at the 12% level for $x\in[0.7,0.9]$ can impact our knowledge of the PDF at $x$ near 1 by more than 20%. This kinematic region is relevant for neutrino experiments, and such precision should be feasible in the near term. In addition, as crosscheck of nuclear theory in this region, exploratory calculations of nuclear PDFs will become available; see Sec. III.4. Such precision is already relevant to neutrino-nucleon scattering at 4–5 GeV. Further, it allows a crosscheck of the nuclear-theory treatment and of the systematic uncertainties of nuclear PDFs.

Neutrino DIS can be important for determining the strange quark and antiquark parton distributions. Currently, no calculation of the Bjorken-$x$ dependence of the strange PDF has been done with lattice QCD, due to numerical limitations, but there are USQCD proposals to investigate $s(x)-\bar{s}(x)$. On the other hand, the nucleon sea flavor asymmetry $\bar{u}(x)-\bar{d}(x)$ has been studied Chen et al. (2018). Unfortunately, the uncertainties in the quasi-PDF approach are currently much larger than those from experimental/phenomenological extraction.

Going to larger nucleon boost momenta $P$ with high statistics is key to reducing several systematic uncertainties in these quasi-PDF and pseudo-PDF approaches Detmold et al. (2019), especially for the antiquark distribution and small Bjorken $x$. However, this poses several computational challenges. First, large momentum translates into large $(aP)^{n}$, and therefore, ensembles with increasingly smaller lattice spacings $a$ are needed. Given the need to keep the spatial size of the box sufficiently large to avoid significant finite-volume effects, which may be enhanced for some nonlocal matrix elements Briceño et al. (2018b), this increases the computational cost. Second, as the momentum becomes larger, the signal-to-noise ratio degrades, even when using methods such as “momentum smearing”Bali et al. (2016), designed to enhance the contribution of the lowest-lying state in correlation functions at nonzero three-momentum, thereby increasing the number of measurements that need to be made. Last, the excited-state contributions themselves become more significant at higher momenta both through the greater number of contributing states arising from the reduced symmetries at nonzero momentum, and through the relative compression of the energy spectrum. This requires either calculations at many source-sink separations, or the use of the variational method with an expanded basis of operators. Thus calculations of the precision that the neutrino program demands requires a significant commitment of computational resources.

III.4 Axial currents in light nuclei

An important challenge for lattice QCD is to extend the calculations of the axial properties of the nucleon to the more complex systems of nuclei. Just as for the nucleon, knowing nuclear matrix elements of axial currents and of quantities relevant in deep-inelastic scattering on nuclei are high priorities in the lattice-QCD community. Over the last decade, first studies of a range of nuclear properties have been performed, and calculations of the requisite axial structure of light nuclei are eminently feasible in a five-year timeframe.

Nuclear effects in neutrino-nucleus scattering are important, extremely nontrivial, and not simply related to those in electron-nucleus scattering. For example, Gamow-Teller transitions in nuclei Buck and Perez (1983); Krofcheck et al. (1985); Chou et al. (1993), which flip spin and isospin of the nucleus, are poorly described in most simple nuclear models, with deviations of as much as $25\%$ from naive expectations based on simply scaling from the single-particle $n\to pe\bar{\nu}$ transition. Once sophisticated nuclear wave functions and many-body axial currents Baroni et al. (2016); Krebs et al. (2017) are used, however, agreement with the data is reached Pastore et al. (2018). At higher energies, state-of-the-art Green-function Monte-Carlo (an exact many-body method) calculations Lovato et al. (2015, 2018) show that neutrino response functions describing scattering on nuclei such as 12C have effects from two-body currents at the 20% level, particularly in the transverse response.

In the last few years, USQCD collaboration members have performed lattice-QCD calculations of $A=2,3$ axial transitions in the forward limit using unphysically heavy quark masses corresponding to $M_{\pi}\sim 800$ MeV Savage et al. (2017); Chang et al. (2018b); Tiburzi et al. (2017). This work, in which the Gamow-Teller contribution to tritium $\beta$-decay and the rate of the $pp\to de\bar{\nu}$ fusion process were extracted, demonstrates that the calculations relevant to constrain neutrino interactions with light nuclei can be performed. While nuclear effects in the axial matrix elements of two- and three-body systems are found to be at the few percent level, the lattice-QCD calculations were able to resolve the relevant effects by isolating the intrinsically two-body contributions. Pursuing these calculations at the physical quark masses and controlling all sources of systematic uncertainties in them will require exascale computing resources. Beyond these forward limit calculations, extensions of this approach to enable calculations of the form factors of light nuclei from nonforward transition matrix elements are underway. Calculations involving multiparticle final states are also necessary but challenging: a theoretical understanding of the simplest inelastic channels is presented in Ref. Briceño and Hansen (2016). For high-energy neutrino-nucleus scattering in the deep-inelastic regime, lattice-QCD calculations of the moments of the relevant parton distributions in nuclei will be useful in constraining nuclear effects in a very different kinematic regime.

While these calculations do not directly address the particular nuclear targets for DUNE and other neutrino-scattering experiments, they are useful in constraining the low-energy constants and meson-exchange currents that enter nuclear-chiral-EFT axial currents Baroni et al. (2016); Krebs et al. (2017). Such matchings are already underway for spectroscopy Barnea et al. (2015); Bansal et al. (2018); Contessi et al. (2017) and electromagnetic interactions Kirscher et al. (2017) at unphysical values of the quark masses, and the machinery necessary to undertake this at the physical quark masses is being developed. As well as studies of currents, lattice-QCD calculations of nuclei and other systems such as three and four neutron systems will provide input into the three-body forces in nuclear EFTs, particularly those aspects of these forces that are challenging to access experimentally.

IV Extremely challenging calculations

Looking further into the future, we can foresee the need for calculations that go far enough beyond the current state of the art that it is hard to know when or whether they will be possible.

IV.1 Electromagnetic and isospin-breaking effects

Going beyond the leading order, calculations of nucleon matrix elements must incorporated the neglected contributions from both QED and from strong isospin breaking (SIB). There are two possible approaches for completing these tasks, which can be generally classified as being perturbative and nonperturbative. Perturbative calculations make use of existing isospin-symmetric and QCD-only lattice ensembles to compute the desired matrix elements. The QED effects are computed with explicit vector current insertions and some scheme for the virtual photon lines. Similarly, the SIB effects are included via scalar current insertions that allow for expansion in quark mass about the isospin-symmetric point. For nucleon matrix elements, perturbative calculations require five-point correlation functions as well as disconnected diagrams and and pose a significant challenge to pursue.

Nonperturbative calculations make use of gauge ensembles that include explicit sea effects for QED and SIB. These are $\text{SU}(3)\times\text{U}(1)$ gauge field ensembles with the up and down quark masses tuned to their physical values. Rather than being restricted to a perturbative expansion of photon fields, these calculations include combined gluon, photon, and quark loops to all orders. While these calculations are likely to be cheaper than the perturbative calculations mentioned above, they also have technical difficulties that must be overcome. The most challenging of these difficulties is likely to come from finite size effects. Since photons mediate a long range force, calculations including QED will be sensitive to the size of the lattice. Many computations on lattice ensembles with different volumes will be necessary to quantify and remove this systematic effect.

IV.2 Axial currents in heavier nuclei

A holy grail for neutrino-nucleus scattering is controlled QCD calculations of the axial form factors, resonance transition form factors, and nuclear PDFs for 40Ar, the target nucleus in DUNE and several other experiments. As yet, exponentially hard challenges must be overcome in order for meaningful lattice-QCD calculations of such heavy nuclei. Both the factorial growth complexity of many-body contractions and the exponentially degradation of the signal-to-noise ratio are currently impede to progress on this front. Since lattice-QCD studies of nuclei are relatively new, it is not unlikely that new algorithms will definitively alter this picture (algorithms involving machine learning and quantum computation Roggero and Carlson (2018) may dramatically improve the situation),555For discussion of these novel approaches, see the companion whitepaper “Status and Future Perspectives for Lattice Gauge Theory Calculations to the Exascale and Beyond” Joó et al. (2019). but at present it is realistic to assume that direct lattice-QCD calculations of argon will not occur in a timeframe relevant for the coming long-baseline experiments.

Perhaps more realistically, significant tests of nuclear EFT frameworks beyond the few-body sector would be enabled by lattice-QCD calculations of the spectrum and axial structure of an intermediate nucleus such as 12C. Aspects of coherent scattering off nuclei will also be addressed by such calculations. While still challenging, a number of groups are investigating ways to perform the relevant contractions and studying improved ways to extract signals from noisy multi-baryon data through optimization methods Detmold and Endres (2014) or improved estimators Beane et al. (2015); Wagman (2017); Wagman and Savage (2017a, b); Detmold et al. (2018). For carbon targets, experimental scattering data exists and comparison of lattice-QCD calculations with this will help understand the systematics of the $A$ dependence of nuclear EFT approaches and assess the reliability of the extrapolations to argon.

For light nuclei, adapting the techniques discussed above to address the Bjorken-$x$ dependence on nuclear PDFs will become possible as computing resources increase. While challenging, and still in the development stage even for the nucleon, these PDF methods will help constrain the flavor and spin dependence of nuclear PDFs that are important for high-energy $\nu A$ scattering.

V computing needs

As we have seen in Secs. III and IV, many topics pertaining to lattice QCD for neutrino physics are still exploratory. In those cases, computing estimates are impossible because progress depends on innovation and flexible computing, rather than an industrial resource. It is feasible and reasonable, however, to estimate the computing needs of the calculations discussed in Sec. II. We do so here, focusing on the example of nucleon form factors.

The methodology for the calculation of the axial and the electromagnetic form factors is well established, and data with control over statistical errors have been generated by several collaborations worldwide. Table 1 list several of these efforts. Unfortunately, only a few include several ensembles with strange sea quarks Ottnad et al. (2018); Gupta et al. (2017); Chang et al. (2018a). Even those calculations should be pushed to smaller lattice spacing and (in one case) smaller up-down quark mass. Furthermore, most calculations obtain mean-squared charge radii ($r_{A}$, $r_{E}$, and $r_{M}$) that are smaller than phenomenological extractions, by about 30%. To diagnose where the difference lies, it is crucial to improve control over systematic uncertainties in lattice-QCD calculations to obtain a definitive result from QCD.

Here we base the cost estimates for achieving a 10–15% result on two ongoing efforts within USQCD. A convenient starting point is the work of PNDME collaboration that has presented extensive results for the axial form factors Gupta et al. (2017) using the Wilson-clover formulation for the valence quarks on ensembles with 2+1+1 sea quarks with the staggered formulation. One way to avoid this “mixed action” approach is to use staggered valence quarks Meyer et al. (2016b). Based on current running on institutional clusters at BNL and Fermilab, we estimate 9M GPU-hours666For QCD codes, these 9M GPU hours correspond to 300M conventional (e.g., Intel Skylake) CPU core-hours. to carry out a calculation on eleven ensembles, five at physical sea-quark mass, and six with $m_{l}={\textstyle\frac{1}{2}}(m_{u}+m_{d})=0.2m_{s}$, with lattice spacing as small as 0.03 fm.

Similarly, a significant subset of USQCD plans on generating a suite of ensembles with Wilson-clover sea quarks. To generate eight such ensembles, with light sea-quark masses corresponding to pion masses of 170 MeV and 270 MeV (four ensembles each), with lattice spacing as small as 0.05 fm. The estimate to finish generating these ensembles is 8M GPU-hours (assuming the GPUs on Summit at ORNL). This chunk of computing will be shared with many other projects, particularly those described in the companion white paper “Hadrons and Nuclei” Detmold et al. (2019). The computation of the needed nucleon correlation functions is estimated to require 15M GPU-hours.

These estimates set the scale for a modern calculation of the simplest quantity needed for neutrino physics. At the same time, comparably demanding work with small nuclei, but not yet physical pion mass, will be needed. Such work is necessary to understand the technical issues facing more realistic calculations and to find better methods and algorithms. Even assuming gains from innovation, it is hard to imagine that nuclear form factors will end up below 10M GPU-hours. The same line of reasoning can be applied to other calculations discussed in Sec. III.

Acknowledgements.

We would like to thank Raúl Briceño, Maxwell Hansen, Richard Hill, Ciaran Hughes, William Marciano, Saori Pastore, Noemi Rocco, and Michael Wagman for useful input. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract No. DE-AC05-06OR23177. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the United States Department of Energy, Office of Science, Office of High Energy Physics.