What makes lattice QCD computationally tractable?

Question

I don't know anything about lattice QCD. Therefore, at first glance it seems to me that lattice QCD should be computationally intractable for all practical purposes.

Let's assume that we only care about the two lightest quark types and their corresponding antiquarks. Then the QCD state should be a joint wavefunction of a certain number of each of the 16 gluons (8 colours times 2 spin states), 12 quarks (2 flavours times 3 colours times 2 spin states), and 12 antiquarks. Except that we also have to sum over all combinations of particle numbers, up to some cutoff where the probability densities become negligibly small.

For example, if there's one of each type of particle, that's 40 particles and each one has 3 coordinates, so our wave function is a function from 120 real numbers to one complex number. Now even if we just divide a volume of space into a $4 \times 4 \times 4$ grid, there are $4^{120}$ possible combinations of inputs to the wave function. So a "solution" for the ground state would have to specify the probability density at $4^{120}$ points. But again, this needs to be summed over all combinations of numbers of particles too, each of which gives an exponent equal to 3 times the total number of particles in the combination.

Anyway, obviously this is not the representation that lattice QCD works with. I am reminded of density functional theory from quantum chemistry, in which the first Hohenberg–Kohn theorem states that all observables of the ground state of the system can be computed from the density functional alone, which is a function mapping 3 spatial coordinates to total electron number density at the corresponding point. It is not necessary to compute the joint wavefunction of all $N$ electrons in the system (which would be a function of $6N$ real numbers, taking into account the two spin states of the electron). Is there some equivalent to this in lattice QCD, e.g., you only need to compute, in each grid cube of the lattice, the number density of each quark flavour, in order to compute e.g. the total energy of the state?

Answer 1

Lattice QCD is used with the path integral formulation, so we don't directly deal with wavefunctions. Consider for now a pure $\mathrm{SU}(3)$ gauge-theory, i.e. no fermions. We now want a (euclidean) action on a lattice $\Lambda = a I_n^4$ for $I_n = \{0, 1, \dots, n-1\}$ that converges to the known Yang-Mills action in the continuum limit $a \to 0$. One possible starting point is the Wilson action \begin{align*} S_\mathrm{Wilson} = \sum_{x\in \Lambda}\sum_{\mu\nu} \frac{1}{g^2}\mathrm{Re}\,\mathrm{tr}\left[1 - P_{\mu\nu}(x)\right] \to \int\mathrm{d}^4x\ \frac{1}{2g^2}\mathrm{tr}\left[F_{\mu\nu}F^{\mu\nu}\right] = S_\mathrm{YM} \end{align*} Here $P_{\mu\nu}(x)$ stands denotes a 'plaquette' which is a square loop of Wilson lines $W$/parallel transporters on the gauge group \begin{align*} P_{\mu\nu}(x) = W(x \to x + a e_\mu) W(x + a e_\mu \to x + a e_\mu + a e_\nu) W(x + a e_\mu + a e_\nu \to x + a e_\nu) W(x + a e_\nu \to x) \end{align*} The Wilson lines are simply the action of finite parallel-transport on the gauge group, i.e. the path ordered exponential of the gauge field, since the gauge field is akin to a connection. This has a nice side effect: in the path integral we typically integrate over the gauge field, i.e. over the whole Lie-algebra, but now we're dealing with finite transformations or elements in the Lie-group which is in our use cases compact. Making the naive integral finite/integrable and we don't have to deal with the whole Faddeev-Poppov construction.

Using this action we can write down the path integral and calculate correlation functions numerically via \begin{align*} \langle A(U)\rangle = \int\mathcal{D}W\ \frac{e^{-S_\mathrm{Wilson}(W)}}{Z} A(W) = \int\left(\prod_{x \in \Lambda}\prod_\mu \mathrm{d}W_\mu(x)\right) \frac{e^{-S_\mathrm{Wilson}(W)}}{Z} A(W) \end{align*} where $W_\mu(x) = W(x \to x + a e_\mu) = W_\mu^a(x) T^a$ can be expanded using the generators of the group. Finally this means that we have a (finite) $D= n^4 \cdot 4 \cdot 8$ dimensional integral. Computing this exactly is as you might guess still unreasonable, but we can use Monte-Carlo methods to sample field configurations according to the distribution $\frac{e^{-S(W)}}{Z}$.

A naive implementation of the procedure above would randomly generate a new Wilson line $W'_\mu(x)$ and see if the action is lower than the action which had $W_\mu(x)$. If the action increased then we might still accept the proposal $W'$ to include finite temperature/quantum corrections. The correlation functions can then be found by averaging over the configurations that were generated this way. This means that per update we only need to generate $n^4$ random $\mathrm{SU}(3)$ matrices and then we need as many iterations of the update process as possible for the best results.

Nb: the Wilson action is not unique, there are infinitely many terms that contain the lattice spacing $a$ that can be included in the action, without changing the continuum limit.