All About the Neutron from Lattice QCD

I describe how simulations of lattice QCD using the path integral formulation provide the two basic quantum mechanical properties of QCD, its ground state in which correlation functions are calculated, and Fock state wavefunctions between which matrix elements of operators are calculated. Both constructs are stochastic, so unfortunately one gets no intuitive picture or even a qualitative understanding of what they look like, nevertheless they contain and display all the subtleties of the quantum field theory. Today, these simulations provide many quantities that are impacting phenomenology and experiments. I illustrate the methods and the steps in the analysis using, as examples, three observables: the isovector charges of the nucleon, the contribution of the quark's intrinsic spin to the nucleon spin, and the pion-nucleon sigma term.


Lattice QCD
The field of lattice QCD (LQCD)-the theoretical rigorous path integral formulation of QCD discretized on a 4D hypercubic grid by Wilson, and used to provide non-perturbative predictions [1]-has come of age.This formalism converts quantum field theories into statistical mechanics systems, for example, the 3+1 dimensional QCD in Minkowski time becomes, after a Wick rotation to Euclidean time, a classical system of gluon and quark fields on a 3+1 dimensional lattice in Euclidean time (see Fig. 1).Numerical simulations of it [2] are providing first principle results with control over all systematic uncertainties for a large number of physical observables that elucidate the standard model and probe physics beyond it.The Flavor Lattice Averaging Group (FLAG) 1 provides a community based evaluation of quantities that are considered robust [3,4]. 2 With improvements in numerical algorithms and increasing computing resources, the errors on these quantities are being reduced steadily, and many more quantities are being added to the list.At the same time, the need for new ideas for a big leap forward is also evident.In this writeup, I will present an idiosyncratic mixture of topics, starting with explaining what simulations of LQCD give us, and then highlight successes and, at the same time, the need for new ideas for sub-percent precision predictions of the properties of nucleons.
!:  "# ! $ :  "# "  % :  "# #  & :  "# $ t A B where g is the gauge coupling and λ a are the Gell-Mann matrices.This lattice formulation preserves gauge invariance of the continuum theory.(Right) Illustration of the path intergral formulation of quantum mechanics of a particle moving between points A and B. Each path i has coordinate x i at time t and is weighted by e iA i , where A i is the action.All possible paths connecting A and B contribute.
I begin with a brief recap of the path integral formulation of a particle moving in time.Figure 1 (right) illustrates some of the paths that contribute to the quantum mechanical amplitude for the particle to go from point A to point B. In fact, all paths one can draw between those two points contribute.Each path n has a weight e iA n , where A n is the action of that path.The value of an observable, say the position x at a prescribed time t is given by the expectation value n x n e iA n / i e iA n where x n is the position at time t.The amplitudes interfere and, typically, the path with the smallest action gives the largest contribution.In the classical limit, this path converges to that predicted by Newton's equations.Using this example, I now motivate how simulations of lattice QCD give us the vacuum configurations, the analogues of the "paths", the correlation functions corresponding to measurements such as x, and the wavefunctions within which matrix elements (ME) of operators (→observables) can be calculated.
In simulations of lattice QCD (LQCD), and of gauge field theories in Euclidean time in general, the analogue of the paths are gauge configurations.Each LQCD configuration C i is a specification of the 12 independent entries in each SU(3) matrix, U x,µ , assigned to each link of the lattice (see Fig. 1).The C i have a weight e −A i , where A i is the Euclidean QCD action calculated on configuration C i .It is a functional of all the U x,µ .Note, the quark fields are formally integrated out as discussed below, leaving only gauge fields as dynamical variables.Configurations are generated using Markov Chain Monte Carlo Methods with importance sampling and the Metropolis accept/reject step [2].Conceptually, this algorithm for the generation of the C i is the same as used in the classical simulations of spin models, however, simulations are computationally expensive because on a 100 4 lattice, there are 4 × 12 × 10 8 independent variables that specify a C i (entries in all the SU(3) link matrices) and evaluating the action A i is expensive.The full set of the C i (called an ensemble) and their associated A i provides us the fully quantum mechanical ground state of the field theory, albeit stochastically since only a finite C i are sampled in practice.The action A is characterized by the input parameters of the simulations: quark masses m i with i ∈ {u, d, s, c} flavors, lattice spacing a (equivalently the gauge coupling β = 6/g 2 via dimensional transmutation) [2], and the lattice volume L 3 × T .
The full set of configurations is the same independent of the action, A(m i , a), i.e., of all input parameters {m i , a}, and depends only on the number of gauge links (4 × L 3 × T ) and the values they can take.There are ∞ 2 of them: infinite number of variables in the limit the volume L 3 × T → ∞, and each variable is continuous valued between {−1, 1}.If even a significant subset of these were needed to calculate observables, precision would not be achieved.

Correlation functions and observables
Expectation values of observables O are obtained from ensemble averages, i Γ α e −A i / i e −A i , of correlation functions Γ α measured on the C i .What saves us from having to consider the ∞ 2 configurations is that the weight e −A i is so very highly peaked about the minimum of A that 10 3 − 10 7 (depending on the observable) importance sampled and statistically independent configurations suffice to yield expectation values with sufficient precision.The location of the peak of the distribution (minimum of A) changes with the input parameters, i.e., with the {m i , a, L, T }.
In practice, data (expectation values of correlation functions) are obtained on many ensembles with different {m i , a, L} so that the limits a → 0, m i to their physical values set by the experimental values of the masses M π , M N , M Ω and M D , and L → ∞ can be taken to obtain physical results.Typical in current simulations, m s and m c , being sufficiently heavy, are already tuned to their physical values before starting production, and only m u,d are varied.In the isospin symmetric limit, m u = m d , it is typical to represent the common light quark mass m ud by the corresponding value of M π .Again, m ud , is tuned before starting production runs.Today, we can perform simulations at M π = 135 MeV, but very often data are also obtained at a number of heavier values of m ud (equivalently M π ), and then extrapolated to M π = 135 MeV using ansatz motived by chiral perturbation theory.These ideas will be illustrated by the calculations/results reviewed later.

021.3
An essential simplification, in fact one that allows simulations of QCD on classical computers in the first place, is that the fermion action for each flavor q, A F = ψDψ = ψ(γ µ (∂ µ + i gA µ ) + m q )ψ → ψ(x + a μ)(γ µ U x,µ )ψ(x) + m q ψ(x)ψ(x) , (1) is bilinear in the quark fields.Here D is the Dirac operator, which on the lattice is a (3 × 4 × L 3 × T ) 2 complex matrix that depends only on the U x,µ .The fermions can, therefore, be integrated out exactly from the path integral but contribute the determinant of the Dirac matrix, DetD f , for each flavor f to the Boltzmann weight, which becomes where A G is the gauge action.They, therefore, continue to impact the A(m i , a, L), and thus the position of the peak of the distribution specifying the ground state and fluctuations about it.Calculating the contribution of the determinant to the Boltzmann factor used to generate the configurations makes simulations expensive but does not pose a formal obstruction.The inverse of D f , a sparse matrix, is the all-points-to-all-points Feynman quark propagator.One column (or row) of D −1 f is the point-to-all-points Feynman propagator S F used to construct correlations functions.
To construct the 2-and 3-point correlation functions of operators composed of quark fields, even though we have integrated them out formally, consider the time-ordered product T of the pion interpolating operator ūγ 5 d and the axial current dγ µ γ 5 u: with the assumption that all three (actually only 2 since momentum is conserved in LQCD) operators have been projected to zero momentum for simplicity, thus leaving only the time index.The notation 〈 • • • 〉 implies ensemble average over the C i .At the same time as integrating out the quarks, one can perform a Wick contraction of the fields in Eq. 2 to get where S F (a column of the inverse of the Dirac matrix D calculated using iterative Krylov solvers) is the Feynman propagator from a point source to all lattice points.Now using the hermiticity property of the Dirac action and its inverse, S F (0, τ) = γ 5 S † F (τ, 0)γ 5 , we get Thus performing the Wick contraction replaces the quark fields in correlation functions in terms of S F = D −1 , which depends only on the gauge links.In short, both in the generation of the configurations and in the calculation of correlation functions, the quark fields are integrated out exactly.The expressions in Eq. 4 for the pion correspond to the quark line diagrams shown in Fig. 2, whose expectation values give the desired non-perturbative correlation functions.
The question I hope you are dying to ask is how does the quark propagator, S F , calculated on a given configuration and combined to form the quark-line diagrams shown in Fig. 2, know anything about the non-perturbative propagation of the pion or any of the thousands of possible states of QCD, the analogous quark-line diagrams for which are obtained by simply changing the interpolating operators.As already explained, the interpolating operators create states with given quantum numbers, for example dγ 5 u or dγ 0 γ 5 u for the pion.These are propagated in time by the transfer matrix.The non-perturbative properties of the propagating pion and its dynamics arise from the coherent addition of those gauge fluctuations on each configuration that correspond to a pion propagating.The miracle of the ensemble average is that only fluctuations corresponding to states with pion quantum numbers survive.In short, the frothing vacuum has all possible fluctuations present, and the ensemble average picks up those that conform to the quantum numbers of the state created by a given interpolating operator.The ensemble average in LQCD simulations creates the stochastic Fock state wavefunction at each time t, indicated by the pink band, and the operator causes transitions between the various "pion" states of the transfer matrix.The extraction of these matrix elements (and thereby pion's axial form factors) are obtained from fits to these 2-and 3-point correlation functions using Eq. ( 5).

Another interesting aside is that
4 is a positive definite quantity.So you may ask-how can averaging over only a small "working ensemble" give a precise unbiased result?The answer lies in the fact that configurations importance sampled according to the Boltzmann weight e −A provide an unbiased approximating to the full ensemble (path integral).Clearly, to improve statistical precision, one needs to enlarge the "working ensemble".Now we come to the last part of the introduction to LQCD-how does one get physics from correlation functions such as those in Eq. 4? For this we invoke the spectral decomposition of Γ 2 π and Γ 3 π , i.e., the insertion of a complete set of "pion" states |π i 〉 at each intermediate time step, and the evolution between time steps given by the transfer matrix.The result is where π is the pion interpolating operator, and the sum over {i, j} is over all the states of the Transfer matrix with the quantum numbers of the pion.Such decompositions of Γ n hold for all interpolating operators and the states they couple to.Simply replace the symbol π by the state of interest.By fitting Γ 2 π versus τ, we can extract the amplitudes, |〈0| π|π i 〉| 2 and the energies E i for all the "pion" states that couple to π.In the limit τ → ∞, only the ground (lowest) state contributes, and for π = dγ 4 γ 5 u, one gets from Γ 2 the pion decay constant , and its mass M π .Thus 2-point functions give us the amplitudes for creating the state and the spectrum of the theory (actually, in discrete time, of the Transfer Matrix).
Next, consider Γ 3 π .It has an additional operator, Âµ , sandwiched between the pion creation and annihilation operators.The propagating tower of "pion" states, |π i 〉, interact with current Âµ at time t, which causes transitions between these states with strength given by the matrix element, 〈π i | Âµ |π j 〉.These ME can be isolated from the fit to Γ 3 π since all the other terms, can, in principle, be determined from the fit to Γ 2 π .It is easy to check that as τ → ∞, the ME within the ground state of the pion is given by the ratio Γ 3 /Γ 2 .
If any two operators in Eq. 4 are projected to ⃗ p = 0 (LQCD conserves momentum), then we get the axial charge of the pion.If the axial current Âµ inserts momentum ⃗ q and one of the pion interpolating operator removes it, we get the axial form factor describing the semileptonic decay of pions with Euclidean momentum transfer Once data for Γ 2 and Γ 3 are collected at a number of values of {m i , a, L}, and fit using Eq. 5 to get data for decay constants, energies E i and matrix elements, their physical valus are obtained by a simultaneous extrapolation: M π → 135 MeV, a → 0, and L → ∞ using physics motivated ansätz.This extrapolation is common to all LQCD calculations as illustrated below.

Renormalization of lattice operators
We can write down a number of equally good lattice operators O latt,n that should give the same results in the continuum limit.At finite a, the results will differ due to their relative normalization and different discretization errors, over and above the known differences coming from the amplitudes, such as 〈π| πi |Ω〉 if different interpolating operators πi are used.Lattice renormalization factors, Z latt O n , relate the different O n at a given a, and their scaling behavior as a → 0. Results with renormalized operators, say Z latt µ,n , should agree in the continuum limit.The experimental results presented by phenomenologists typically use a scheme such as M S and a convenient scale such as 2 GeV above which perturbation theory is considered reliable.To translate the lattice result to the M S scheme at, say, 2 GeV is a two step process.First one calculates the lattice factors Z i O in some scheme (currently two popular ones are the regularization independent [symmetric] momentum schemes labeled RI-MOM and RI-sMOM [5,6]), and a second calculation that relates them to M S that is typically done in the continuum using perturbation theory, as is the factor for running in the continuum to a specified scale, say, 2 GeV.
For the calculations described here, the renormalization factors are well-determined.Many other operators, such as the CP-violating Weinberg and quark chromo EDM operators of dimension 6 and 5, respectively, have divergent mixing with lower dimension operators.Cnstructing finite renormalized versions to use in simulations is very non-trivial [7].In fact, for these two operators, it is still an open problem.

Nucleon correlation functions
The quark line diagrams for the nucleon 2-and 3-point functions are shown in Fig. 3. Formally, the mechanics of the lattice calculation is very similar to that for the pion, however there are two very important differences: • The signal to noise ratio falls exponentially as ∼ e −(E N −1.5M π )τ in all nucleon correlation functions, whereas the pion has no degradation.Typical data show that for Γ 2 N a good signal extends to about 2 fm and for Γ 3  N to about 1.5 fm with O(10 6 ) measurements [8].
• For a number of matrix elements, excited states contributions (ESC) from towers of multi-hadron excited states, N π, N ππ, . . .labeled by their relative momentum, are enhanced and still large at 1.5 fm [8].Their energies begin at about 1200 MeV, much below radial excitations.Fully removing these ESC in fits to the spectral decomposition of Γ 2 and Γ 3 remains a challenge for many observables.
To determine various quantities, we use appropriate probes.Changing the operator to a scalar, Ŝ = du or tensor, T = dσ µν u, gives us nucleon's scalar and tensor charges that are also probed in precision measurements of neutron decay distributions [9].One link operators give us the momentum fraction, helicity and transversity moments [10].And the list continues.
Having laid out, hopefully, an intuitive introduction to simulations of lattice QCD, I now discuss three calculations in order of increasing complexity.

Isovector charges of the nucleon
The iso-vector axial, scalar, tensor charges of the nucleon, , and g u−d T , probed in the N → P decay, are extracted from Γ 3 (⃗ p = 0), i.e., from the forward matrix element with Dirac matrix X O = γ µ γ 5 , , σ µν specifying the insertion of the axial, scalar and tensor operators at zero momentum transfer.For these iso-vector charges, only the connected quark line diagram (middle panel in Fig 3) contributes in the isospin symmetric limit, and results for the proton and the neutron are the same.The data in Fig. 4 for g u−d A from a {a = 0.071fm, M π = 170MeV} ensemble (see Ref. [8]) illustrate what the presence of ESC does and our goal is to understand and reliably remove them.The data (same in the two panels) display the following features of the ESC: • The variation of the data with t and τ is the signature of ESC.In the limit τ − t and τ → ∞, the data (Γ 3 /Γ 2 = g u−d A ) should be flat in t and lie on top of each other, i.e., independent of τ and t, particularly near t − τ/2, i.e., away from the source and sink.
• The data should be symmetric about t − τ/2 because Γ 3 is.The statistical quality of the data for τ = 19 (=1.35 fm) is already borderline in this respect.
• The convergence of the data with τ for fixed t is monotonic and from below.This shows that ESC cause g u−d A to be underestimated.
• ESC is removed by fitting data at the 3 largest values of τ using the spectral decomposition given in Eq. 5 truncated at 3 states (ground plus two excited states).The value for the ground state matrix element, given by the fit, is shown by the grey band.
• The data in the two panels are the same.The fits differ in the energy, E 1 , of the first excited state used.In the left panel it is the output of the fit to Γ 2 while in the right, the energy of the lowest N π state with relative momentum (0, 0, 1) is input using a narrow prior.The motivation for including this N π state is χPTit contributes at one loop.
• The value of the ground state ME given by the two fits is different but the augmented χ 2 /d o f are comparable.This highlights a serious problem: current data are not at sufficiently large τ nor precise enough for fits to discriminate between different choices of E 1 .• The two first excited-state energies, E 1 , selected are, physics wise, reasonable options: the first is given by the fit to Γ 2 and lies close to the N(1440), while in the second fit we input, N (0, 0, 1)π(0, 0, −1).Furthermore, N (0, 1, 1)π(0, −1, −1) and the rest of the tower also contributes.In fact, all states with the same quantum numbers contribute!What we do not know, a priori, are the amplitudes, and the size of the contribution of each possible excited state to the ME.In short, the statistical precision of the current data allow fits with three states, however, these fits show that there are large regions in E 1 , E 2 and A i that give similar χ 2 /d o f but significantly different τ → ∞ values.

021.7
Bottom line: Until the data are good enough to distinguish between fits with different number or combinations of plausible excited states, and lacking a theoretical reason for a particular choice, the difference between the extrapolated values with different possible excited states can be regarded as an estimate of the systematic uncertainty due to ESC.A factor of 10 increase in statistics will give data for τ = {21, 23} with precision similar to τ = {17, 19} data shown in Fig. 4.Then, I believe, fits with different E 1 will start to be discriminated by χ 2 /d o f .
Our data suggest that the uncertainty due to including or not the N (0, 0, 1)π(0, 0, −1) state could be a ∼ 5% effect in g u−d A , but is much smaller in g u−d S or g u−d T [8].At this point in time, controlling ESC remains the key outstanding systematic for achieving sub-percent precision in the prediction of the isovector properties of nucleons.
The next step in the analysis, illustrated using the g u−d A data, is the chiral-continuum-finitevolume (CCFV) extrapolation using a simultaneous fit in {M π , a, M π L} shown in Fig. 5.The methodology is the same as described in [8], except we now have data on 13 ensembles with different {M π , a, M π L}. (Note, m ud is specified by M π .)For the three iso-vector charges, the extrapolation ansatz, keeping lowest order corrections in each of the 3 variables, is [8]: In such ansätz, the order of the correction terms one needs to include depends, in general, on the level of improvement of the lattice action and the observable.Each panel in Fig. 5 shows the result of this fit versus a single variable with the other two set to their physical values.For example, in the panel versus a, the pink band shows the result with M π set to 135 MeV and 021.8 M π L → ∞.The data, however, have not been shifted in these two variables, which is why they do not lie within in the pink band.
Such simultaneous CCFV fits are routine for getting physical values of observables.Features in these fits that illustrate the size of the three systematics are • For the 2+1-flavor clover-Wilson action we have used, the value of g u−d A decreases by about 10% as a → 0, i.e., the slope with respect to a is positive.

• The value of g u−d
A increases by about 10% over the range of M π as M π → 135 MeV.
• The dependence on a and M π is largly independent of each other, with opposite slopes.
• There are no significant finite volume corrections seen for M π L > 4.This welcome feature has been observed in all calculations involving single nucleon correlation functions.
Phenomenologically, the most interesting of the isovector charges is the axial charge, g u−d A , which has been extracted from experiments with high precision, g u−d A /g V = 1.2754(13) [11].The precision and robustness of lattice results have increased over the last decade, but my conclusion, in light of possible unresolved ESC of multihadron states such as the N π that can cumulatively be as large as ≈ 5%, is we need to better quantify and remove the ESC contributions before we can claim sub-percent precision.

Contribution of the spin of the quarks to the nucleon spin
Using Ji's gauge invariant decomposition 1/2 = q=u,d,s,c ∆q/2 + L q + J g [12], where L q is the quark orbital and J g the gluon total angular momentum, the contribution of the intrinsic spin of a quark with flavor q, ∆q/2, to the proton spin is given by a relation very similar to Eq. 6: with ∆q = g q A and u and ū are the quark spinors.Because the operator is diagonal in flavor, there is now an additional Wick contraction in which the operator forms a closed loop as illustrated in the right panel in Fig. 3.This is called a "disconnected diagram".The full contribution to g q A is the sum of the connected (middle) and disconnected (right) quark line diagrams.The calculation of disconnected diagrams introduces a new layer of computational cost.The straightforward solution to calculate the momentum projected quark loops with operator insertion on all time slices is to calculate the all-to-all quark propagator D −1 .This is not practical as it is a (12•10 8 )×(12•10 8 ) complex matrix for a 100 4 lattice.The solution is to construct  [3,4] obtained by different collaborations ( PNDME '22 [14] ETM '19 [15], PNDME '18 [16] Mainz '18 [17], χQCD 18 [18], JLQCD '18 [19], ETM '17 [20], χQCD 15A [21], Engelhardt '12 [22]).The points with filled squares meet the FLAG criteria for inclusion in the FLAG average.The contribution to the nucleon spin from quarks with flavor q is ∆q/2 = g q A /2. a stochastic estimate.This approach works well, is bias-free but introduces additional statistical uncertainty in the sum due to the stochastic estimation of the disconnected quark loop diagram whereas the calculation of the connected quark-line diagrams is exact up to matrix inversion precision for S F .This uncertainty in the measurement of the sum on each configuration gets convoluted with that due to gauge fluctuations in the ensemble average.Methods such as deflation and bias-corrected truncated solver methods have allowed the reduction in errors in the disconnected contributions to be of the same size as in the connected, however, since their central value is smaller they contribute a larger fraction to the overall error.
The steps in the analysis to get g u A , g d A , g s A are the same as for the isovector charges described in Sec. 5. Our final results (PNDME 2018 and 2022 (preliminary)) are shown in Fig. 6 along with those from other collaborations in the FLAG format [3,4].In the FLAG review process, results that pass the criteria for control over discretization, finite lattice volume, renormalization and ESC systematic uncertainties, and obtained sufficiently close to physical pion mass (or extrapolated to M π = 135 MeV), are then averaged and the overall error estimated with appropriate consideration given to possible correlations between results.These FLAG averages [3,4] are the community consensus value.

The pion-nucleon sigma term
The pion-nucleon σ-term, σ πN , is a fundamental parameter of QCD that quantifies the amount of the nucleon mass generated by up (u) and down (d) quarks having non-zero mass.In the "direct" method discussed here, it is, in the isospin symmetric limit m ud = (m u + m d )/2, given by where the extraction of m ud and g u+d S are done separately.The scalar charge g q S is determined from the forward matrix element of the scalar density qq between the nucleon state: (SM): it determines the coupling of the nucleon to any scalar mediator with quark content of the coupling given by qq, for example in direct-detection searches for dark matter (DM) scattering off nuclei via a scalar mediator (similarly g q A give the spin dependent and g q T the tensor couplings); in lepton flavor violation in µ → e conversion in nuclei; and in electric dipole moments.The calculations of g q S and g q T are similar to that discussed above for g q A . Figure 7 shows data for g u S + g d S from a physical pion-mass ensemble [23].Again the two fits have very similar χ 2 /d o f but the one with E 1 of N π as the excited state gives almost 50% larger value.Our NNLO χ P T analysis given in Ref. [23] shows that there are two significant ESC, one from N π and the other from N ππ, and each contribute about 10 MeV to σ πN .These two states are almost degenerate in our lattice calculation, so they effectively contribute as one in the fit to Γ 3 , i.e., their amplitudes in the spectral decomposition add as the exponential factors are very similar.The right panel in Fig. 7 shows pictorially why the N π state makes a large disconnected contribution: the scalar current has a large coupling to the quark loop, and the configuration shown in the quark-line diagram favors a N π intermediate state.
The two analyses with different values of E 1 led to an interesting conundrum.The standard analysis (with E 1 ∼ 1450 MeV) gave σ πN ≈ 40 MeV consistent with previous lattice analyses [3,4], whereas the analysis with E 1 = E N π ≈ 1230 MeV gave σ πN ≈ 60 MeV, which is consistent with the dispersive analysis starting with the N π scattering data [23].Our preferred solution is the latter based on the χ P T analysis.In that case the tension between LQCD and phenomenological estimates is resolved.
Clearly, the 50% difference between the two analyses with similar χ 2 /d o f calls for additional LQCD calculations to be done to confirm this exciting result.The key point for future calculations of σ πN , using either the direct method of calculating the charges g u,d S as defined in Eq. 9 or using the Feynman-Hellmann relation, σ πN = m q ∂ M N /∂ m q where M N is the nucleon mass, is that they have to be done close to M π = 135 as only there the N π state becomes much lighter than N (1440) and the ESC are very different and manifest.Our data with M π > 200 MeV do not give significantly different results between the two kinds of fits.Thus, extrapolation from heavier M π ensembles will miss this physics in both methods.

Conclusions
Simulations of LQCD provide ensembles of importance sampled configurations whose distribution according to the Boltzmann factor DetDe −A G = e −A G +LndetD constitutes the nonperturbative ground state of QCD.This construction is exact (bias-free) but stochastic.For a given lattice action, ensembles are characterized by six input parameters {m i , a, M π L} with m i ∈ {m u , m d , m s , m c }. Correlation functions of any time-ordered string of operators are given by quark-line diagrams obtained using the Wick contraction.Properties of QCD (spectrum, matrix elements, EoS, etc) are extracted from expectation values, i.e., ensemble averages, of these correlation functions.The full excursion (possible "paths") of the quark propagators (and values of link parameters for non-local and gluon operators) in the quark line diagrams over 3-space but at a fixed intermediate time generate the full Fock space wavefunction, a linear combination of all states with the same quantum numbers as the interpolating operator.The propagation of each of these states in Euclidean time τ is damped as e −E n τ , allowing the ground state with energy E 0 to be isolated in the limit τ → ∞.Once again, this stochastic description of the wavefunction provides no intuition or visualization but allows the calculation of fully quantum mechanical matrix elements of any operator within this state.
Lattice QCD is called a "black box" because we cannot visualize or represent the vacuum fluctuations or the wavefuctions that get created at intermediate times in the correlation functions Γ n .Nevertheless, the results obtained are rigorous, display all the subtleties of QCD and confirm that this quantum field theory describes the quantum dynamics of quarks and gluons.
Three kinds of observables have been used to illustrate how LQCD calculations are done and data analyzed.Looking ahead, precision calculations of nucleon correlation functions need to overcome two challenges: the exponential degradation of the signal and how to remove all/most excited state contamination in the wavefunctions to get matrix elements within the nucleon ground state from correlations functions calculated with finite source-sink separation τ.The ESC in current data can be large as illustrated by the pion-nucleon sigma term.
The future is exciting -the FLAG reports provide a growing testimony to LQCD having matured and results having an impact on phenomenology and experiments [3,4].Methodology and algorithms for many calculations are robust, however, brute force approach to achieving sub-percent precision in nucleon correlation functions and ME by just increasing the statistics is unlikely to succeed in the next few years.It is, therefore, time for innovation and an exciting challenge to the next generation-develop new methods and algorithms to reduce systematics and increase statistics efficiently.

Figure 1 :
Figure 1: (Left) Discretization of QCD on a hypercubic lattice with quark fields placed on sites and the gluon fields A µ (x) on directed gauge links via SU3 matrices U x,µ = e ia gA a µ λ a

FockFigure 2 :
Figure 2: Illustration of quark-line diagrams for 2-point (left) and 3-point functions for the pion (middle and right).The gluon lines are just for illustration and to remind the reader that all orders of gluon exchanges are implicit in these diagrams.(Right) The axial current, shown by , is inserted at intermediate Euclidean time t and with momentum ⃗q.The ensemble average in LQCD simulations creates the stochastic Fock state wavefunction at each time t, indicated by the pink band, and the operator causes transitions between the various "pion" states of the transfer matrix.The extraction of these matrix elements (and thereby pion's axial form factors) are obtained from fits to these 2-and 3-point correlation functions using Eq.(5).

χ 2 /χ 2 /Figure 4 :
Figure 4: for the ratio Γ 3 A /Γ 2 with the insertion of the axial current with ⃗ p = 0 at time t are shown with different colors for 4 values of the source-sink separation τ. Results for the ground state (τ → ∞) ME (= g A ), with two different plausible mass gaps, a∆M 1 , in a 3-state truncation of the spectral decomposition of Γ 3 A , are shown by the grey band, and for different values of τ by lines of the same color as the data.The two estimate of g A are different but the fits are not distinguished by the χ 2 /d o f .

Figure 5 :
Figure 5: Result of the simultaneous chiral-continuum-finite-volume fit, Eq. 7, to g u−d A data is shown by the pink band and plotted versus the lattice spacing a (left), pion mass M 2 π (middle), and M π L, the lattice size in units of M π (right).The grey band in the middle panel shows a simple chiral fit, i.e., with c 2 = c 4 = 0 in Eq. 7.

Figure 7 :
Figure 7: Fits to get the scalar charges g u S + g d S using for the mass gap a∆M 1 the value obtained from Γ 2 (left) versus the noninteracting energy of the N π state (middle).The right panel shows the disconnected diagram, whose contribution is large.