{'fraction_non_alphanumeric': 0.035142105400265716, 'fraction_numerical': 0.01325969729335452, 'mean_word_length': 4.615564657694558, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Twenty-first Century Lattice Gauge Theory: Results from the QCD Lagrangian

Quantum chromodynamics (QCD) reduces the strong interactions, in all their variety, to a simple nonabelian gauge theory. It clearly and elegantly explains hadrons at short distances, which has led to its universal acceptance. Since its advent, however, many of its long-distance, emergent properties have been believed to be true, without having been demonstrated to be true. This paper reviews a variety of results in this regime that have been established with lattice gauge theory, directly from the QCD Lagrangian. This body of work sheds light on the origin of hadron masses, its interplay with dynamical symmetry breaking, as well as on other intriguing features such as the phase structure of QCD. In addition, nonperturbative QCD is quantitatively important to many aspects of particle physics (especially the quark flavor sector), nuclear physics, and astrophysics. This review also surveys some of the most interesting connections to those subjects.


Introduction
Quantum chromodynamics (QCD) is the modern theory of the strong nuclear force. It is part of the Standard Model of elementary particles and the underpinning of terrestrial and astronomical nuclear physics The conception of QCD is rightly hailed as a triumph of reductionism, melding the quark model, the idea of color, and the parton model into a dynamical quantum field theory. At the same time, the scope of QCD is rich in emergent phenomena. Symmetries emerge in idealized limits: C, P , and T are exact when the total "vacuum angle"θ = 0; chiral symmetries emerge when two or more quark masses vanish (1,2); and heavy-quark symmetries are revealed as one or more quark masses go to infinity (3,4). More remarkable still are the phenomena that emerge at a dynamically generated energy scale Λ QCD , the "typical scale of QCD." Much of what is known about QCD in this nonperturbative regime has long been based on belief. Evidence from high-energy scattering fostered the opinion that QCD explains the strong interactions and, therefore, the belief that QCD exhibits certain properties; otherwise, it would not be consistent with lower-energy observations. These emergent phenomena-such as chiral symmetry breaking, the generation of hadron masses that are much larger than the quark masses, and the thermodynamic phase structure-are the most profound phenomena of gauge theories. The primary aim of this review is to survey how lattice QCD has enabled us to replace beliefs with knowledge. To do so, we cover results that are interesting in their own right, influential in a wider arena, qualitatively noteworthy, and/or quantitatively impressive.
The rest of this article is organized as follows. Section 2 introduces the QCD Lagrangian and discusses how, in a general setting, to fix its free parameters. Section 3 gives a short summary of lattice QCD methodology. Sections 4 and 5 discuss hadron masses and their connection to chiral symmetry. An output of these calculations are the quark masses and the gauge coupling, which are discussed in Section 6, along with some timely results pertaining to flavor physics. Section 7 presents some interesting properties of nucleons. Section 8 discusses the phase structure of QCD. Section 9 offers some perspective. The appendix identifies resources for readers who wish to start research in numerical lattice gauge theory.

Quantum Chromodynamics
The (renormalized) Lagrangian of QCD has "1 + n f + 1" free parameters (where n f is the number of quark flavors): where F µν is the gluon's field strength, / D = γ µ (∂ µ + A µ ), and ψ f denotes the quark field of flavor f . The first parameter is the gauge coupling g 2 ; the next n f parameters are the quark masses m f ; and the last parameter,θ, multiplies an interaction that violates CP symmetry. Experiments have demonstrated the existence of n f = 6 quarks. At energies below the top, bottom, and charm thresholds, however, it is convenient and customary to absorb the short-distance effects of these quarks into a shift of g 2 and then take QCD with n f = 5, 4, or 3.
The coupling g 2 diminishes gradually with increasing energy, stemming from virtual processes of gluons and the n f active quarks; this "running" is known as asymptotic freedom (5,6).
In the Standard Model, quark masses arise from the matrix of Yukawa couplings to the Higgs field, y. The matrix y can be brought into a form with real eigenvalues, y f = √ 2m f /v, and an overall phase, arg det y. (Here, v = 246 GeV is the Higgs field's vacuum expectation value.) In this context, the coupling multiplying ε µνρσ tr[F µν F ρσ ] is altered:θ = θ − arg det y, where θ is allowed in QCD as soon as CP violation is admitted. Only the differenceθ is observable.
Before one can state that a mathematical theory describes or explains the natural world, one must fix its free parameters with the corresponding number of measurements, in this case 1 + n f + 1. Because the color of quarks and gluons is confined, the free parameters of QCD must be connected to the properties of QCD's eigenstates, which are the color-singlet hadrons. From this perspective, the parameters of QCD may be fixed as follows. The electric-dipole moment of the neutron is too small to measure, which leads to a boundθ < 10 −11 . Such delicate cancellation of θ and arg det y is a mystery, known as the strong CP problem (7). For QCD calculations it means simply that we can setθ = 0 with no significant consequences. The rest of the parameters are tuned to reproduce 1 + n f specific hadronic properties. Because the gauge coupling runs, the physical interpretation of g 2 is predicated on the energy at which it reaches a fiducial value, say, g 2 = 1. But the mathematics is, strictly speaking, dimensionless, so the energy at which g 2 = 1 obtains a physical meaning as a multiple (or fraction) of some standard mass, such as the proton mass. In practice, it is wiser to choose the standard mass to be insensitive to the quark masses. Finally, the n f quark masses are best related to n f hadron masses that depend sensitively on them; for example, the kaon mass is used to tune the strange-quark mass, because M 2 K ∝ m s . With lattice gauge theory (8), one has a tool to relate the QCD Lagrangian directly to such hadronic properties and, thereby fix the parameters. Hadronic properties, by the way, always fix the parameters of QCD. The top quark mass, for example, is measured at the Tevatron via the four-momenta of hadrons in jets.

Numerical Lattice QCD
Lattice gauge theory (8) was invented in an attempt to understand asymptotic freedom without gauge-fixing and ghosts (9). The key innovation is to formulate nonabelian gauge invariance on a space-time lattice. Then the functional integrals defining QCD correlation functions are well defined, because DU , Dψ, and Dψ are products of a countable number of individual differentials. Here S = d 4 x L QCD is the action, • is just about any gaugeinvariant product of fields, and Z ensures that 1 = 1. This formulation is formally equivalent to classical statistical mechanics, which allows theorists to apply a larger tool kit to quantum field theory. For example, Wilson (8) used a strong-coupling expansion to lowest order in 1/g 2 to demonstrate confinement. The results presented below have been obtained by integrating expressions of the form (2) on big computers with Monte Carlo methods. Lattice gauge theory defines QCD mathematically and, thus, in principle provides an algorithm for computing anything. Nevertheless, the computer imposes practical constraints. To compute anything within a human lifetime, the integrals are defined at imaginary time, t = −ix 4 , which turns Feynman's phase factor into the damped exponential of Equation 2. A computer, obviously, has finite memory and processing power, so the spatial volume and time extent of the lattice must be finite.
These limitations do not impair the computation of many important classes of quantities. The imaginary time imposes no problem whatsoever for static quantities. The finite volume introduces errors in one-particle states that are exponentially suppressed and, hence, a minor source of uncertainty. In twoparticle states, the finite-volume effects are stronger, but the volume dependence yields information such as scattering lengths. Similarly, effects of the finite time extent are exponentially suppressed, except in thermodynamics, where it becomes a tool. Finally, the continuum limit must be taken as part of the renormalization procedure (10,11).
From Equation 2, it is straightforward to derive some simple results for correlation functions. The two-point function is where π is a composite field of definite quantum numbers (e.g., of the pion), and the sum ranges over all radial excitations. For a large enough time separation x 4 , a fit to an exponential yields the lowest-lying m π 1 and | 0|π|π 1 |. By using a larger set of operators, one can extend this method to compute excited-state properties. For a transition with no hadrons in the final state, as in leptonic decays, one can simply replace π(x 4 ) with a current J, J(x 4 )π † (0) = n 0|Ĵ |π n π n |π † |0 e −mπ n x 4 (4) in which the only new information is 0|Ĵ |π n , yielding for large x 4 the decay matrix element of the lowest-lying state. For a transition with one hadron in the final state, such as one from a B meson to a pion, one needs a three-point function, in which the new information is π n |Ĵ|B m , yielding for large x 4 , y 4 the matrix element between the lowest-lying states. We compute matrix elements for flavorchanging processes, dark-matter detection, and nucleon structure in this way. Equations (3)-(5) are derived by inserting complete sets of eigenstates of the QCD Hamiltonian. The only assumption is that these eigenstates are hadrons. Thus, every successful fit of these formulae for hadronic correlators provides a posteriori incremental evidence that hadrons are indeed the eigenstates of QCD.
In all cases of interest, the fermion action is of the formψMψ, where the spacetime matrix M is a discretization of the Dirac operator (plus quark mass). Then the fermionic integration in Equation 2 can be carried out by hand: Here, the fermionic integration replaces ψ iψj in • with [M −1 ] ij to yield • ′ . Importance sampling, which is crucial, is feasible only if det M exp (−S gauge ) is positive.  (14,15) U In most cases, a notable exception being the case of nonzero baryon chemical potential, this condition holds. The determinant det M represents virtual quark-antiquark pairs, also known as sea quarks. The matrix inverse M −1 is the propagator of a valence quark moving through a stew of gluons A and sea quarks. Several quark propagators are sewn together to form hadronic correlation functions, which via Equations 3-5 yield masses and transition matrix elements. The sea det M poses the largest, and the propagator M −1 the second largest, computational challenge. The numerical algorithms become even more demanding as the quark mass is reduced. Lattice QCD data with unphysically heavy up and down quarks can be extrapolated to the physical limit guided by chiral perturbation theory (12,13). This step removes the cloud of unphysically massive pions and replaces it with a better (and improvable) approximation to the physical pion cloud.
Because det M and M −1 require so much computing power, several formulations of lattice fermions are used. As one might anticipate, the computationally fastest and theoretically cleanest methods are not the same. (If one formulation were both fastest and cleanest, no one would use anything else.) Each formulation can be characterized by the amount of flavor symmetry retained by the lattice (Table 1). Staggered fermions are computationally the fastest, but the flavor group comes in a semi-direct product with the symmetry group of the hypercube, SW 4 , and the total number of species in the continuum limit, for n f fields, is 4n f . This fermion doubling is not a problem for the propagators M −1 stag . For the sea, however, one must take (16) [det M stag ] 1/4 and appeal to numerical and perturbative evidence that the rooting yields a local field theory in the continuum limit (19).
Because of the expense of sea quarks, many lattice QCD calculations have been carried out with two or fewer (light) flavors. The error entailed in omitting the strange-quark sea is difficult to estimate, so this review considers mostly results with 2+1 flavors in the sea. Here 2+1 denotes the strange sea and two more flavors, for up and down, taken as light as possible. Such simulations made a breakthrough early in this century (20). Now the first results with the charmed quark sea, "2+1+1," are becoming available.

Hadron Spectrum
We compute the masses of hadrons not only for a quantitative comparison of QCD with nature but also to learn how gauge theories generate mass. As Section 6 shows below, hadron masses are much larger than the sum of the masses of the underlying quarks. The positive binding energy stems from the confining properties of the gluon field and from the kinetic energy of the quarks.
Let us begin with the energy in a flux tube between a static quark-antiquark pair, as a function of their separation r. The lowest energy level of the flux tube is the potential energy V (r), and the excitations of the flux tube are also informative. The states are labeled Σ ± g,u , Π g,u , ∆ g,u , . . ., according to the eigenvalues of gluonic angular momentum along r and of CP [in the subscript g (u) stands for (un)gerade, which is German for even (odd)]. The Σ states also carry a superscript ± denoting the change of sign (or not) of the wave function upon reflection in the plane containing r; otherwise such reflections relate degenerate pairs. Figure 1 shows the lowest-lying levels in the SU(3) gauge theory without light quarks (21). At short distances, the level spacing and ordering is consistent with asymptotic freedom: V (r), for example, is Coulombic up to logarithmic corrections. As r increases, the spacing changes, and at a separation of around 2 fm, the level ordering rearranges to that of a string. The level spacing does not become fully string-like until larger separations (22). The behavior of the excitations is instructive, because the lowest level, V (r), becomes consistent with a string at a relatively short distance around 1 2 fm (23). A vivid picture of the flux tube has it narrowing as r increases, owing to the attraction between gluons, but the details suggest that the flux tube looks more like a sausage than a string.
One can imagine connecting the ends of the sausage to obtain non-qq states known as glueballs. Such mesons have no counterpart in the quark model, and lattice gauge theory provides the best (theoretical) evidence that these states do indeed exist. Glueball masses with 2+1 flavors of sea quarks show little change (24) from earlier calculations with no sea quarks (25). In particular, the masses remain consistent with the idea, motivated by lattice QCD, that the f J (1710) is the is the heavy-quark potential, exhibiting Coulombic behavior at short distances and linearly confining behavior at large distances. The higher excitations also exhibit the level ordering of electrodynamics at short distances but the level ordering of a string at large distances. The colors pink, orange, green, and violet stand for decreasing lattice spacing. Data are from Reference 21.
Lattice QCD has been used to verify the mass spectrum of quark-model hadrons within a few percent. Figure 2 shows four broad efforts on the spectrum of the isopsin-1 light mesons and the isospin-1 2 and -3 2 baryons (27,28,29,30,31). All these simulations include 2 + 1 flavors of sea quarks, and the error bars in References 27, 28, 30 reflect thorough analyses of the systematic uncertainties. A satisfying feature of Figure 2 is that the results do not depend in a systematic way on the fermion formulation chosen for the quarks. Even the latest results for the difficult η-η ′ splitting are encouraging (32,33,34). Figure 2 includes predictions for mesons with quark contentbc (38, 36, 39). The prediction for the pseudoscalar B c has been (subsequently) confirmed by experiment (40,41), whereas the prediction for the vector B * c awaits confirmation. These predictions build on successful calculations of the bb and cc spectra (37,42,43,44,45), which reproduce the experimental results well.
The most striking aspect of the spectrum is how well it agrees with nature. The nucleons provide almost all the mass in everyday objects, and their masses have been verified within 3.5%. Their mass mostly comes, via m = E/c 2 , from the kinetic energy of the quarks and the energy stored in the sausage-like flux tube(s) holding the quarks together.

Chiral Symmetry Breaking
A striking feature of the hadron spectrum in Figure 2 is that the pion has a small mass, around 135 MeV, whereas the other hadrons have masses more than five times larger. To understand the origin of the difference, Nambu (1) applied lessons from superconductivity, noting (four years before quarks were proposed) that a spontaneously broken axial symmetry would constrain the pion's mass to vanish, with small explicit symmetry breaking allowing it to be nonzero. If the up and down quarks can be neglected, the QCD Lagrangian acquires an SU L (2) × SU R (2) symmetry, thereby explaining the origin of Nambu's axial symmetry. The consequences of spontaneous symmetry breaking were studied further by Goldstone (46), which led to the formula (47), when applied to QCD with massless up and down quarks. The flavor-singlet vacuum expectation value ψ ψ is known as the chiral condensate. If either factor on the left-hand side of Equation 7 is nonzero, then the other must vanish.
From the earliest days of QCD, most physicists were confident that this theory explained the richness of the strong interactions. Because both QCD and Nambu's picture of the pion were considered correct, it was believed that QCD must generate a chiral condensate. Until recently, however, no ab initio calculation of ψ ψ tested Equation 7. Lattice QCD has now filled this gap (48,49,50): Here, the first uncertainty is statistical and the second is a combination of systematics, and the quark masses have been adjusted to Nambu's idealization, m u = m d → 0, m s physical (50). In summary, QCD's symmetries and dynamics have now been demonstrated to account for the hierarchy between the pion and the other hadron masses.

Standard Model Parameters
The Standard Model (with nonzero neutrino masses and mixing angles) has 28 free parameters: • Standard electroweak symmetry breaking: v = 246 GeV, Lattice QCD is essential or important in determining the values of eleven parameters (the first under gauge couplings and all but m t under quark sector). Lattice field theory (without QCD) is also useful for shedding light on the Higgs self-coupling λ (51) and the top Yukawa coupling y t = √ 2m t /v (52). As in Equation 8, this section reports several results in the conventional MS scheme at a renormalization point µ, such as 2 GeV or the Z boson mass M Z . MS is the most common renormalization scheme in quantum field theory, because it is the technically simplest one in perturbative calculations. Readers unfamiliar with the MS scheme should focus on the consistency and on the relative uncertainty of the results, rather than the details of the scheme's definition.  These results are remarkable for at least two reasons. First, the up and down masses are very small, approximately four and nine times the electron mass, respectively. Quark masses arise from interactions with the Higgs field (or something like it). Thus, this sector is not the origin of much of the mass of everyday objects. Second, m u , although very small, is nevertheless significantly far from zero. This outcome is interesting because if m u = 0, then the additional symmetry of the Lagrangian would renderθ unphysical, obviating the strong CP problem (7). (A subtlety could arise from a nonperturbative additive correction to m u , but it is probably too small to alter this conclusion.) The heavy charmed quark and bottom quark masses can be determined along the same lines, but the most accurate results come from computing quarkonium correlation functions and taking their continuum limit (58). These functions give spacelike information on the same function measured in the timelike region in e + e − → cc or e + e − → bb. Perturbation theory to order α 3 s (59) can then be used to extract the heavy-quark masses and α s (57). This approach yields the results in the fifth column of Table 2, which are in nearly perfect agreement with the corresponding determinations from e + e − collisions (60).
Lattice QCD provides excellent ways to determine the gauge coupling α s = g 2 /4π. In lattice gauge theory, the bare coupling g 2 0 is an input. Alas, for many lattice gauge actions, perturbation theory in g 2 0 converges poorly (61), obstructing a perturbative conversion to the MS or other such schemes. Two other strategies have been adopted to circumvent this obstacle. One is to compute a shortdistance lattice quantity-such as a Wilson loop-and reexpress perturbation theory for the Monte Carlo results in a way that eliminates g 2 0 in favor of a renormalized g 2 (61,62). The other is to compute a short-distance quantity with a continuum limit, and then apply continuum perturbation theory. The  Table 3 and compared with an average of determinations from high-energy scattering and decays (70). There is excellent consistency among the results, not only with different discretizations of the determinant for sea quarks but also when the charmed sea is included (69). An important source of uncertainty is the truncation of perturbation theory, including strategies for matching to the MS scheme, and running to scale M Z . In the example of the small Wilson loops, an independent analysis of the data from Reference 66 has been carried out; this analysis found α s (M Z ) = 0.1192 ± 0.011 (71), which should be compared with the second line of Table 3.
The agreement between the lattice QCD and perturbative QCD results for α s , m c , and m b is especially compelling because QCD is the union of the quark model of hadrons and the parton model of high-energy scattering. This consistency is evidence that the QCD of hadrons and the QCD of partons are the same.

Flavor Physics
As mentioned in Section 2, the quark masses arise from the electroweak interactions. In a basis where the mass matrix is diagonal, the W boso couples to all combinations of (u, c, t, . . .) ⊗ (d, s, b, . . .) T quarks. Along with the SU(2) gauge coupling, the udW vertex carries a factor V ud , and similarly for all other combinations. As a change of basis (from the weak interaction basis to the mass basis), the Cabibbo-Kobayashi-Maskawa (CKM) (72,73) quark-mixing matrix V is unitary. After considering global symmetries of the gauge interactions, one sees that the CKM matrix has fewer parameters than a generic unitary matrix does. For three generations, three mixing angles and one CP -violating phase remain, and a convenient choice consists of |V us |, |V ub |, |V cb |, and arg V * ub . Lattice QCD calculations play a key role in flavor physics. The phase and, except for |V tb |, all magnitudes of the CKM matrix can be accessed via processes for which the corresponding lattice-QCD calculations are under good control. (Other processes that do not need lattice QCD also provide information on the CKM matrix.) A representative set of calculations and their utility is illustrated by mixing (86), the full suite of lattice QCD calculations overdetermines the CKM matrix and, thus, tests for consistency. The semileptonic D decays are considered crosschecks. Taking |V cd | and |V cs | from CKM unitarity (which is very precise), one finds that lattice QCD calculations of the kinematic distributions (83) and the normalization of the rate (84,85) agree with results from several experiments. Non-SM particles could spoil this picture, which is why it is interesting to test it in detail. With a fourth generation of quarks and leptons, the 3 × 3 submatrix generically would not be unitary. If other particles, such as supersymmetric partners of the known particles, change quark flavor, then the SM relation between a flavor-changing process and V is spoiled. The off-diagonal elements are small-|V us | ∼ 0.2, |V cb | ∼ 4 × 10 −2 , and |V ub | ∼ 3 × 10 −3 -so it is not out of the question for (widely anticipated) TeV-scale particles to make detectable contributions.
During the first decade of the twenty-first century, all simple lattice-QCD calculations that are pertinent to the CKM matrix were carried out with 2+1 sea quarks. In most cases, more than one collaboration has published results, and, in many cases, the literature covers more than one fermion formulation for the quarks. The calculations most directly connected to determining the CKM parameters have been combined-with an eye to correlations in the errors-in Reference 87. (Updates are available at http://latticeaverages.org/).
Despite the broad agreement between flavor-physics measurements and the Standard Model, some tension appears in global fits to the four CKM parameters (88). Some mild discrepancies also arise in a few isolated processes, and here we consider two leptonic decays in which lattice QCD plays a key role.
Let us begin by noting that the semileptonic and leptonic determinations of |V us |, which rely on the matrix elements of References 74,75 and 90,91,92, respectively, are completely compatible (89). The vector component of the W boson mediates the former, and the axial current mediates the latter. Because these determinations are compatible, nothing other than the SM W boson with its V − A coupling is needed to account for these decays.
The history of results on the D s → µ + ν and D s → τ + ν decays suggests caution, however. In 2008, the measured branching fractions exceeded the SM predictions by nearly 4σ (Figure 3). This discrepancy relies on the decay constant f Ds from lattice QCD (100,90). Through the use of the same methods as for the π and K decay constants (90,101), the D s decay constant can be computed to 1-2%. At the time, some experimenters asserted that either the lattice QCD calculations were wrong or that New Physics mediated the decay. For example, the excesses of D s → ℓ + ν could be due to leptoquarks (102); a few-percent amplitude could be constructively interfering. (Leptoquarks are hypothetical particles with lepton and quark quantum numbers.) As more measurements came in, however, the discrepancy softened, and it is now only 1.6σ. Although the lattice QCD average of the D s decay constant has increased by 8 MeV, the experimental average has decreased by 18 MeV. The early measurements fluctuated upward; perhaps the same holds for B → τ ν.   7 Nucleon Matrix Elements, Dark Matter, and the LHC Two of the most compelling questions facing particle physics are the origin of electroweak symmetry breaking and the composition of dark matter. The experiments developed to address them rely on the familiar proton and neutron. The Large Hadron Collider (LHC) collides pp, and the Tevatron pp, and detectors buried deep underground hope to observe weakly interacting massive particles (WIMPs) scatter of the protons and neutrons in nuclei. To interpret the results of these experiments, it is helpful to calculate certain matrix elements of the nucleon (106). Let us start with WIMP-nucleon scattering. Each quark's contribution to the cross section is proportional to a matrix element known as the sigma term, where m q is the quark mass, and M N is the nucleon mass. All quarks q in the nucleon can contribute, including virtual quarks such as s and c. The partial derivative should be taken with the other n f QCD parameters held fixed. Usually the light quarks are combined into the isospin-singlet (Beware of factors of 2 in the definitions of these and similar quantities, sometimes denoted Σ, in the literature.) For WIMP detection, σ u+d and σ s are especially important, because virtual pairs of heavier quarks are far too uncommon inside nucleons. Until recent lattice QCD calculations became available, the light quark sigma term was extracted from πN scattering, with the help of chiral perturbation theory (χPT) (107, 108). As shown in the first two rows of the σ u+d column of Table 4, the extraction depends more on assumptions than on the experimental data. To estimate σ s , one uses information from the baryon octet masses (109). Unfortunately, this information must be subtracted from σ u+d , which renders σ s rather unstable. There is a pressing need to improve both matrix elements (110).
Equation 10 suggests two methods to compute σ q in lattice QCD: either from a three-point function, as in Equation 5, or through study of the mass dependence of the nucleon mass M N . The latter is known as the Feynman-Hellmann theorem, and here one can either reweight the Monte Carlo ensemble to take   (119), and ETM (120). The last has 2+1+1 flavors of sea quarks, the others have 2+1 flavors. The black squares show two fits (121,122) to experimental data; other recent fits of this kind fall between these two (122).
the derivative locally or study the chiral extrapolation to obtain a global handle on the derivative. Table 4 compiles several recent results. From a quantitative perspective, it seems that the results are still settling down, although they tend to favor lower values of σ u+d . More strikingly (and self-consistently with low σ u+d ), the results for σ s substantially contradict the high values used in WIMP phenomenology. Even if these results are not yet as mature as those reported in Sections 4, 5, and 6, they seem to give a more stable picture than the non-lattice estimates. Further work and promising new methods (117) should clarify the sigma terms soon.
In pp or pp collisions, the essential long-distance ingredient in computing cross sections are the moments of the parton distribution functions. These moments are given by matrix elements of local operators, similar to theqq in Equation 10 but with different Dirac structures, such as γ µ and γ µ γ 5 , and derivatives to pull out higher powers of the momentum fraction. Figure 4 shows some recent lattice QCD results for the nonsinglet average momentum fraction x u−d as a function of simulation m 2 π (118, 119, 120), compared with two phenomenological results (121,122). The latter are obtained by fitting experimental data, which exist over a large but limited range of x, to reasonable parameterizations. In principle, the lattice QCD moments add extra information, but the status of the chiral extrapolation may preclude this step at this time, even though some functional forms lead to results (118,119) that agree with the fits to experiment. In this regard, earlier work with 2 flavors of sea quarks yielded confusing results. In a few years, the low moments of quark densities from 2+1-and 2+1+1-flavor simulations should become good enough to incorporate into the traditional fits of experimental data. For collider phenomenology, however, the real challenge for lattice QCD is to compute similar moments of the gluon density, which are less well constrained by low-energy experiments.

QCD Thermodynamics
The previous sections consider isolated hadrons at zero temperature. Soon after the Big Bang, however, the universe was much hotter than it is now. In neutron stars, for example, the baryon density is much higher than in normal nuclear matter. These phenomena have motivated the study of the thermodynamics of QCD. Even within lattice gauge theory, thermodynamics is a vast subject (123,124), so this review touches only on some of the more fascinating aspects.
Thermodynamics starts with thermal averages in the canonical ensemble where T is the temperature, and the traces Tr are over the Hilbert space of the QCD HamiltonianĤ. In fact, the average on the left-hand side of Equation 12 is precisely that of Equation 2; the time extent N 4 specifies the temperature T = (N 4 a) −1 . The eigenstates ofĤ, which are single hadrons and multiparticle states composed of hadrons, do not change with T , but as T increases, the vacuum no longer dominates the way it does in Equations 3-5, and multihadron states begin to play a role in the thermal average. The simplest observables are quantities suach energy, pressure, entropy density, and order parameters sensitive to symmetry breaking. The thermal state can either restore a spontaneously broken symmetry of the vacuum or be a state of broken symmetry itself. Of course, the (approximate) symmetry of the Hamil-  tonian remains intact. Figure 5 shows order parameters for deconfinement and chiral symmetry restoration, as the temperature increases from normal hadronic matter to a phase known as the quark-gluon plasma. Both order parameters change dramatically for a temperature around 145-170 MeV (125,126), but neither, especially deconfinement, exhibits the sharp change characteristic of a phase transition. Studying a whole suite of thermodynamic observables confirms that the transition is a smooth crossover (127,128). This result came as a surprise, and the next two paragraphs explain why. The crossover means that as the early universe cooled, hot matter gradually became more and more like a gas of distinct hadrons. With a first-order phase transition, on the other hand, bubbles of the hadronic phase would have formed inside the quark-gluon plasma. Without a real phase transition, the quark-gluon plasma is not necessarily a fluid of quasi-free quarks and gluons. The eigenstates in Equation 12 remain color singlets, but a thermal medium can be qualitatively different. First, thermal fluctuations encompass states with many overlapping hadrons, so color can propagate from one hadron to the next, as if deconfined. Second, the thermal average applies nearly equal Boltzmann weights to states of both parities, so chiral symmetry can be restored in the thermal average, even though the vacuum breaks it.
The nature of the QCD phase transition is influenced by the physical values of the up, down, and strange quark masses. For vanishing quark masses, the transition would be first order, but as the masses are increased, the strength of the transition diminishes. As depicted in Figure 6a, the physical quark masses (Table 2) are just large enough to render the transition a crossover. If the light quark masses-crucially m s -were around half their physical size, the universe would have cooled through a first-order transition. Before lattice QCD established these results, the conventional wisdom was that the quark masses are somewhat larger than shown in Table 2, yet small enough to remain in the first-order basin of massless quarks.  At nonzero baryon density (chemical potential µ = 0), the fermion determinant becomes complex, which is an obstacle to importance sampling. This restricts lattice QCD calculations to small µ. It is thought that the transition becomes first order for µ ∼ few 100 MeV (Figure 6b) (129), but the matter is not yet settled (130).

Summary and Outlook
The topics discussed above demonstrate that we have learned a great deal in this century about QCD from lattice gauge theory. The twenty-first century is still young, and the prospects for learning more are bright.
Although the mass spectrum of the lowest-lying hadrons has been well verified, it will be interesting to extend the calculations to excited states (131,132) and even to small nuclei such as the deuteron (133) or the H dibaryon (134,135). Beyond QCD, one may wonder whether nature uses gauge theories to generate quark, lepton, and weak boson masses (136).
Most of the calculations related to flavor physics are entering an industrial phase, where the objective is higher and higher precision. An exception is the measurement of direct CP violation in the kaon system. This calculation requires a two-pion final state, and although the formalism for handling this state has long been available (137), only now have K → ππ amplitudes become feasible (138,139). These finite-volume techniques are related to methods (140) for computing scattering lengths (141,142), which have many applications in hadronic physics.

Appendix: Tools
Research in lattice QCD requires computer time and software. Through several efforts around the world, these needs pose lower obstacles than ever before. Several groups have made documented software available so that new programs can be modeled after existing code, rather than being built from scratch. Furthermore, many groups make ensembles of lattice gauge fields available, principally via the International Lattice Data Grid (ILDG) (http://www.usqcd.org/ildg/). In exchange for a suitable citation of an article describing the content of the ensembles, anyone can use these simulation data for his or her own physics analyses. In many cases, even more ensembles are available from collaborations with newer ensembles under generous terms: Most of these collaborations have some core physics analyses but are happy if the expensive simulation data are mined for more results.
Publicly available software can be obtained from the USQCD Collaboration (http://www.usqcd.org/usqcd-software/). Newcomers should start with one of the applications packages: chroma, cps, milc, or FermiQCD. A useful tutorial on this software has been put together by Joó (144; contains slides only).
Two kinds of computing are important to lattice gauge theory, capability and capacity. One needs access to computers with the greatest capability-those able to run large-memory jobs with huge appetite in CPU time-to generate the ensembles of lattice gauge fields. On these ensembles an analysis consists of a huge number of small-to-medium computing demands; this step requires computers with high capacity. Many university groups have access to computers of sufficiently high capacity to analyze the publicly available ensembles.

Acknowledgments
The author thanks Jimmy Juge, Julius Kuti, & Colin Morningstar for supplying the data plotted in Figure 1 Literature Cited