Bootstrap for Lattice Yang-Mills theory

We study the $SU(\infty)$ lattice Yang-Mills theory at the dimensions $D=2,3,4$ via the numerical bootstrap method. It combines the Makeenko-Migdal loop equations, with a cut-off $L_{\mathrm{max}}$ on the maximal length of loops, and positivity conditions on certain matrices of Wilson loops. Our algorithm is inspired by the pioneering paper of P.Anderson and M.Kruczenski but it is significantly more efficient, as it takes into account the symmetries of the lattice theory and uses the relaxation procedure in line with our previous work on matrix bootstrap. We thus obtain rigorous upper and lower bounds on the plaquette average at various couplings and dimensions. For $D=4$, the lower bound data appear to be close to the MC data in the strong coupling phase and the upper bound data in the weak coupling phase reproduce well the 3-loop perturbation theory. Our results suggest that this bootstrap approach can provide a tangible alternative to the, so far uncontested, Monte Carlo approach.


INTRODUCTION
The SU (N c ) Yang-Mills (YM) lattice gauge theory (LGT) is a fundamental ingredient of modern particle physics.Its most illustrious applications are the Standard Model and, in particular, the Quantum Chromodynamics.Nowadays, most of the non-perturbative computations in Yang-Mills theory are done by Monte Carlo (MC) simulations for the lattice formulation of YM theory.Combined with the perturbation theory (PT) [1][2][3][4][5] and RG tools, MC methods have had a huge success, especially in the recent couple of decades, due to the development of supercomputers.It allowed us to compute with a reasonable precision certain masses of hadrons and the S-matrix elements in QCD, reproducing the experimental data [6,7].However, the absence of any systematic non-perturbative "analytic" alternative to MC is, practically and intellectually, somewhat uncomfortable.Moreover, the MC method has its inherent limitations: statistical errors, finite lattice size, high numerical cost of inclusion of dynamical quarks, difficulties in treating finite baryon density, and the real-time dynamics.
An interesting alternative for the study of the LGT is provided by Makeenko-Migdal loop equations (LE) [8,9] for Wilson loop averages (WA).An early attempt at numerical study of LE in the large N c , 't Hooft limit was proposed in [10][11][12][13], in the form of minimization of effective action in the loop space.A more recent brave attempt to bootstrap the LE, combining them with certain positivity conditions [14,15] revived hopes of a more analytic approach.Slightly later, a similar bootstrap method was proposed in [16] for the multi-matrix models.In our work [17] we significantly improved matrix bootstrap by introducing a "relaxation" procedure and applied it to an "analytically unsolvable" large N 2-matrix model, with remarkable efficiency and precision, noticeably exceeding those of MC for the same model [18].

FIG. 1: Our bootstrap results for upper and lower bounds on plaquette average in 4D
LGT (1): for L max = 8 (yellow domain) for L max = 12 (orange curves) and L max = 16 (blue curves).Red circles represent the MC data for SU (10) LGT (with 5 purple squares for SU (12)).Dashed upper and lower lines represent the 3-loop PT (13) and strong coupling expansion (14), resp.These developments have been considerably inspired by the success of the bootstrap approach to CFTs [19,20] and S-matrices in massive QFTs [21][22][23].
Unlike MC where the result is given up to statistical error bars, the bootstrap methods provides rigorous inequalities giving upper and lower bounds on computed physical quantities.These bounds can only improve when increasing the number of bootstrapped variables and constraints on them.
Here we develop a powerful numerical bootstrap algorithm for solving the LEs [8,9] in the lattice YM theory at N c → ∞ and demonstrate it on the computation of plaquette average u P = 1 Nc ⟨tr U P ⟩.Its main ingredients are i) positivity of correlation matrix of WAs; ii) our FIG.2: The plot ∆u P ≡ u boot P − u PT P which might capture the non-perturbative values of the gluon condensate ⟨tr (F µν F µν )⟩.relaxation procedure of [17]; iii) positivity of reflection matrices due to lattice symmetries.iv) symmetry reduction of the positivity conditions.In the supplementary material we worked out a simple example for our general method together with some data points.
Our data (obtained on a single workstation), already for the modest length cutoff L max = 16, look quite encouraging: as seen on Fig. 1, for D = 4, our lower(upper) bounds are quite close to the MC data [14,15,24,25] in the strong(weak) coupling phase, at least far enough from the phase transition.The upper bound is remarkably close to the 3-loop PT.
On Fig. 2 we plot the difference of the bootstrap upper bound and the 3-loop PT: ∆u P ≡ u boot P − u PT P as function of 1/λ.It might capture the non-perturbative effect for the gluon condensate ⟨tr (F µν F µν )⟩ -in principal a measurable observable [26][27][28].The graph is very smooth and it is positive even slightly beyond the phase transition point λ c ≃ 2.9.[29]

YANG-MILLS LOOP EQUATIONS AT LARGE N
We study the LGT with the Wilson action [30]: where U P is the product of four unitary link variables around the plaquette P and we sum over all plaquettes P , including both orientations.The main quantities of interest in 't Hooft limit N c → ∞ are the WAs: The matrix product goes over the link variables belonging to the lattice loop C. W [C] are subject to LEs, i.e., the Schwinger-Dyson equations expressing measure invariance w.r.t.group shifts U l → U l (1+iϵ).Schematically

LE reads:
ν⊥µ where the LHS represents the loop operator acting on the link l µ by replacing it with the loop around plaquette −−→ δC ν lµ or ←−− δC ν lµ (depending on the orientation, as shown in the first line of Fig 3).The LHS sum goes around all 2(D−1) µν-plaquettes orthogonal to the direction µ.The RHS sum goes over all appearances of the the same lattice link l in the loop C. The RHS product corresponds to splitting of the contour C → C ll ′ •C l ′ l , as explained in the 2nd and 3rd lines of Fig. 3. Finally, ϵ ll ′ = ±1 for links l and l ′ with opposite or colinear orientation, respectively.For more details on LEs, see [14,31,32].

Back-track loop equations
To get the full list of the LEs, we also consider the "back-track" LEs.They correspond to doing variations on the links at the end of "back-track" paths originating from the vertices of Wilson loop.These "back-tracks" are equivalent to inserting the identity, but their Schwinger-Dyson variation, can give independent LEs.Fig 4 shows an example of "back-track" LE: The LEs close on single trace WAs (2) (due to the large N c factorization of color traces), which means that FIG.4: An example of nonlinear "back-track" LE. they behave as classical quantities on the loop space.At L max = 16, we have about 40, 000 LEs in 3D and about 100, 000 LEs in 4D, and the back-track LEs constitute more than 80% of all these LEs.Around three-quarters of them are linearly independent and only a small minority are non-linear.Finding the solution to these equations is our primary task for any other progress in studying the physical quantities of planar QCD or its 1/N c corrections.The problem is very complicated since it is formulated in extremely complex loop space.We will try to solve it using the bootstrap approach.

BOOTSTRAP ALGORITHM Positivity Constraints
In general, our positivity conditions come from the positivity of possible inner products on the vector space or a subspace of the operators, i.e.
where O = α i O i is an operator with arbitrary coefficients α i , and O i are basis vectors of the operators.
One of possible adjoint operators O † comes from taking the Hermitian conjugation [14,16,17].For a Wilson path, the Hermitian conjugation corresponds to reversing the path.By taking a linear combination of all Wilson paths 0 → x (between the points 0 and x), with arbitrary coefficients, we can get non-trivial positivity conditions from their inner product.For example, we have only two paths 0 → (1, 1), at and the positivity condition reads: This gives u 2 P ⩽ 1, obvious from unitarity.We call the positivity matrices arising from the Hermitian conjugation the correlation matrices.
Apart from Hermitian conjugation, we have additional reflection positivity conditions where adjoined operators O † come from reflection symmetries [33].For LGT, there are three types of reflections w.r.t.different planes: site, For computations in this work, we consider the full positivity constraint 0 → x for any possible x when L max ⩽ 12.But for L max = 16, we consider only the paths 0 → 0 for various positivity matrices since: 1.When constructing correlation matrices, all the positivity conditions on the open Wilson paths 0 → x are already contained in 0 → 0 correlation matrix for higher lengths (due to back-trackings).
2. At L max = 16, we observe empirically that the 0 → 0 constraints are computationally the most efficient.One important reason for that is that the positive matrices corresponding to 0 → 0 are numerically more tractable w.r.t.symmetry reduction that we will discuss below.

Convex relaxation
In general, LEs are non-linear.Here we use the relaxation method to replace all the non-linear LE with linear ones where in the r.h.s.we replace the products of WAs with new variables Q ij , subject to extra constraints [17]: Here W = {W 1 , W 2 , W 3 , . . .} denotes the column vector of all inequivalent WAs.Notice that the relaxation matrix has rank=1 precisely when The LEs combined with the convex relaxation and positivity conditions constitute the constraints of semidefinite programming (SDP).To get rigorous bounds on WAs, we can maximize or minimize u P .

Dimension Hermitian
Conjugation

Reduction by Symmetry Group
In conformal and S-matrix bootstraps, the positivity conditions are well-known to be organized in different spin channels [19,36] and different irreducible representations of global symmetry [37][38][39][40].In parallel to that observation, we can greatly reduce the positivity matrices via the global lattice symmetries.
Formally, if we have an invariant group G preserving the inner product then the positivity condition on the matrix M defined in Eq 5 can be re-arranged into a block-diagonal form corresponding to the irreps of G.This well-studied procedure is known under the name "Invariant Semidefinite Programming".Here we refer to a statement [41] directly related to our current problem: If the vector space of the paths can be decomposed as a direct sum of irreducible representations Rep k of the invariant group with multiplicity m k : then the positivity condition of the inner-product matrix is equivalent to the collection of positivity conditions on the matrices corresponding to each Rep k , with matrix dimension m k × m k .
For the correlation matrix 0 → 0, the invariant group For the reflection positivity matrices 0 → 0 the invariant groups are subgroups of B d ×Z 2 , leaving the reflection plane invariant.These invariant subgroups are summarized in Table I.
Implementing this symmetry reduction is similar to projecting the physical state w.r.t spin and parity in conformal or S matrix bootstrap.Practically we do the following steps: 1. Find a specific realization of every irrep of the invariant group using GAP software [42].
2. Use the algorithm initiated in [43] to find an equivalent real representation (if the irrep by GAP is complex).

3.
To decompose into such irreps, we use the projector to Rep k [44] : Here r αβ is a matrix element of a real representation identified at the step 2, and α, β = 1, 2, . . ., dim(Rep k ).Taking α = 1, P k = p 11,k gives us a projector to Rep k .

Selection of multiplets of Wilson paths
The Wilson paths form different multiplets of the invariant group.Within each multiplet, the symmetry group permutes different Wilson paths.When constructing the positivity matrices, some multiplets are more important than others.We kept only the most important multiplets.More precisely, several WAs are not related to other WAs by the LEs, such as the 4×4 square Wilson loop at L max = 16.We believe that such WAs and the open Wilson paths out of which they are constructed are relatively unimportant.
As an example, take the correlation matrix for the paths 0 → 0 at 3D and L max = 16.It has a huge size 6505 × 6505.After the symmetry reduction and selection of the multiplets, the positivity of the correlation matrix reduces to postivity of 20 smaller matrices, each corresponding to its irrep, with sizes: So the SDP gets greatly simplified.

DISCUSSION OF RESULTS
Here we present the results of computation of plaquette average u P (λ) for various bare couplings λ in LGT in 2D, 3D, and 4D.On Fig 1 we compare our bootstrap data in 4D with MC for SU (10) [15] and SU (12) [24] LGT, assuming that, with our accuracy, N c = 10, 12 are close enough to N c = ∞.We also compare it with the known 3-loop N c = ∞ PT result [46]: as well as with the SC expansion valid in the SC phase, beyond the 1st order phase transition point λ c ≃ 2.9 [47]: ). ( 14) FIG.6: u P for 2D (upper) and 3D (lower) LGT: the upper and lower bounds from our bootstrap at L max = 8 (yellow region), L max = 12 (orange curves) and L max = 16 (blue curves).The 3D and L max = 12 result without reflection positivity (gray curve) is much less constraining.The line of red circles represents the MC data for SU (10) LGT.The dashed black curve in 2D plot is the exact solution (15).The dashed black curve in 3D plot is the 3-loop PT result [45].
The bootstrap bounds on u P given for L max = 8, 12, 16 on Fig. 1 are quickly improving with the increase of cutoff.The physically most interesting WC phase is much better described by the upper bound.Moreover, we see that the upper bound nicely reproduces the 3-loop PT (13) for a large range of coupling, even beyond the phase transition point where PT is, strictly speaking, not valid.However, comparing these results to MC data, we see that it is not yet so good at capturing the departure from the PT in the interval (2.4, 2.8) where the MC data of [24] (given by black squares on Fig. 1) were used to compute the masses and the string tension.We expect a significant improvement for this range in our data if we reach L max = 20 or even 24.However, this will certainly demand much bigger computational resources.
Finally, we briefly discuss 2D and 3D cases.For 2D LGT the plaquette average can be computed exactly [48,49] (as well as any loop average, see [50]): This example was important for both checking our algorithm and for observing how fast our bootstrap data approach the exact result when increasing L max .The results are presented on Fig. 6.For physically interesting and challenging case of 3D LGT, we compare on Fig. 6 our bootstrap bounds at L max = 8, 12, 16 with the MC data [15] as well as with the known 3-loop PT [45] and SC [47] results.We observe a reasonably fast approach of bootstrap bounds to the MC data when increasing L max , but they are not as close to PT as in 4D case.We employ in LEs the WAs up to the maximal length L max , so it can be considered as our IR cutoff.The physical scale l ph (set by inverse mass or square root of string tension) should, ideally, satisfy 1 ≪ l ph a L ≪ L max , where a L is the lattice spacing.In this paper, we have, for the best of our data, L max = 16 which suggests that the window for the scale of measurable physical quantities in lattice units should be roughly 2 ≲ l ph a L ≲ 6 (compare it to Table 8 of [24] where the IR cutoff is set by the size of space-time torus, typically in the range 10 to 16, and the typical physical length is set by string tension 3 ≲ l ph a L ≲ 6).We conclude that, even though the currently achieved values of L max in our bootstrap approach may be not sufficient to match the precision and scope of the MC experiments, especially for the confinement sensitive physics (glueball masses, string tension, etc.), our results give good hopes on a considerable improvement when augmenting L max .Moreover, for sufficiently small couplings our upper bound data are already at least as good as MC.

PROSPECTS
The bootstrap procedure proposed here has a clear perspective for improving our results and advancing towards the computation of interesting observables.In particular, by choosing objective functions other than u P , one can hope to get a better estimate for all involved physical quantities.For our current implementation at L max = 16, every data point takes ∽ 20 hours of CPU time for 4D, and only half an hour for 3D [51].First, we want to increase the cutoff to L max = 20 and even to L max = 24.This will certainly need supercomputer power.From our current results, we expect a quick narrowing of our bounds to the accuracy comparable to, or even better than MC (without its toll of statistical and systematic errors).Furthermore, since in the 't Hooft limit we don't have internal fermion loops, we can try to find the quark condensate and hadron masses by simply summing up the WAs with the spinorial factors for the relevant one-and two-point functions.The 1/N c corrections (which might be small enough even for the physical N c = 3 case) seem to be not insurmountable tasks since they are subject to linear LEs [9], with coefficients given by the solution of LE (3).One of such problems is the computation of glueball masses from the connected correlator of two small Wilson loops.One can also try to bootstrap directly the N c = 3 YM theory, where the absence of large factorization could be compensated by multiple functional relations between WAs, absent for N c = ∞ case.
Supplementary Material for "Bootstrap for Lattice Yang-Mills theory" Vladimir Kazakov, Zechuan Zheng A worked-out example This section gives a step-by-step solution of L max = 8 bootstrap bound for 2-dimensional lattice Yang-Mills theory in the large N c limit.For simplicity, here, we will only take into account positivity conditions from Wilson lines starting from the origin and ending at the origin.In this situation, up to rotations and reversing the line, we have only one multiplet under symmetries: Taking the inner product between the different Wilson lines in this multiplet, we can construct the positivity matrices corresponding to Hermitian conjugation (correlation matrices) and space reflection (reflection matrices).For the correlation matrix, we have: The row above the matrix is the Wilson lines for the inner products, and the column left to the matrix is their Hermitian conjugation.The I denotes the length 0 trivial Wilson line.The inner product is defined by joining the path and its Hermitian conjugation.This matrix contains six independent Wilson loops: (S2) We denote them as: respectively.(Here, the W 1 is the u P in the main text.) We also note to the reader that although at length 8, there are Wilson loops like the 2 × 2 Wilson loop doesn't appear in the correlation matrix.We will also see that this is true for the loop equations and the reflection matrices, i.e., this will be an SDP with six variables.
There are two reflection matrices corresponding to site reflection and diagonal reflection.In terms of W variables, they read: The reflection matrix corresponds to link reflection is trivial here.Although not evident, at the present bootstrap cutoff, the positivity of reflection matrices won't improve the constraint.So we won't further discuss these matrices (including their symmetry reduction) and do the optimization only under the positivity of the correlation matrix.
Finally, we do the symmetry reduction on the correlation matrix.This process will make the matrix factorize to a series of smaller matrices corresponding to each irreducible representation of the global symmetry group.We stress that this step is purely group theoretical, it doesn't change the positivity constraint itself, but it makes our algorithm orders of magnitude faster.
In our current dimension, as shown in the main text, the global symmetry is actually B 2 × Z 2 .The B 2 group is the cubic discrete version of O(2), and isomorphic to the Dihedral group D 4 .The Z 2 group corresponds to the symmetry of reversing the Wilson loop.The Z 2 group is the Charge Conjugation symmetry on the lattice system, ensuring that all the Wilson loop averages are real numbers.
The Wilson lines above the matrix Eq S1 actually form a multiplet for the global symmetry group.We can decompose this multiplet into each irreducible component.And from the discussion in the main text, the Wilson lines in different irrep are orthogonal in the corresponding inner product.
For the group B 2 × Z 2 , there are ten irreducible representations: The reader could verify that these new positivity conditions are equivalent to the Eq S1.Still, as SDP conditions, they are much more (several orders of magnitude for higher order bootstrap problems) efficient.By the loop equation Eq S5 and the factorized positivities Eq S9, we can maximize or minimize W 1 to get the upper and lower bounds.For λ = 2, we have: We stress that all the loop equations are linear at this lowest order of bootstrap.We must implement the convex relaxation described in the main text to deal with the quadratic terms for higher order bootstrap.

Selected four-dimensional data
In this section we collect some selected data from the 4D bootstrap plot.

-FIG. 3 :
FIG. 3: Schematic representation of LEs: The first line shows the variation of a link of Wilson loop in the LHS of Eq.(3).The 2nd and 3rd lines show the splitting of the contour along the varied line into two sub-contours, for two different orientations of coinciding links in the RHS of Eq.(3).

FIG. 5 :
FIG. 5: Examples of three reflection symmetries on the lattice allowing new positivity conditions on WAs.
where B d is the Hyperoctahedral group in d spacetime dimensions.It acts on a Wilson path by corresponding rotations and reflections on the spacetime lattice.Z 2 is the group action reversing the path.

TABLE I :
Invariant groups of correlation and reflection matrices.
We note that for the two irreps correspond to two E, there are other two vectors, but they will give the same positivity matrices.They are fully redundant.Taking the inner product between the vectors in the same irrep, we get the following positivity condition: S7)The first index is the standard notation for irreducible representation of Dihedral group D 4 ⋍ B 2 , and the second one is the parity of Charge Conjugation.Using the algorithm stated in the main text, we can project the multiplet Table II compares the bootstrap upper bound on u P with the perturbation and MC data in the weak-coupling phase at several specific couplings.

TABLE II :
Plaquette average at various 't Hooft couplings λ: bootstrap, perturbation theory and MC results.