Loop Equations and bootstrap methods in the lattice

Pure gauge theories can be formulated in terms of Wilson Loops correlators by means of the loop equation. In the large-N limit this equation closes in the expectation value of single loops. In particular, using the lattice as a regulator, it becomes a well defined equation for a discrete set of loops. In this paper we study different numerical approaches to solving this equation. Previous ideas gave good results in the strong coupling region. Here we propose an alternative method based on the observation that certain matrices $\hat{\rho}$ of Wilson loop expectation values are positive definite. They also have unit trace ($\hat{\rho}\succeq 0, \mbox{tr} \hat{\rho}=1$), in fact they can be defined as density matrices in the space of open loops after tracing over color indices and can be used to define an entropy associated with the loss of information due to such trace $S_{WL}=-\mbox{tr}[ \hat{\rho}\ln \hat{\rho}]$. The condition that such matrices are positive definite allows us to study the weak coupling region which is relevant for the continuum limit. In the exactly solvable case of two dimensions this approach gives very good results by considering just a few loops. In four dimensions it gives good results in the weak coupling region and therefore is complementary to the strong coupling expansion. We compare the results with standard Monte Carlo simulations.


Introduction
Gauge theories are of fundamental importance for our understanding of Nature but many of their properties are still mysterious, for example in pure gauge theories, the phenomenon of confinement is still not fully understood. Even further, the AdS/CFT correspondence [1] has shown that gauge theories in the strongly coupled regime can be described equally well in terms of string theory in a higher dimensional space. That means that certain gauge theories contain quantum gravity, emergent space-time and strings as bound states. Fundamental to this understanding is the relation of gauge theories and string theory in the limit of a large number of colors as envisioned by 't Hooft [2]. In such relation, exemplified by AdS/CFT, the string theory description is completely in terms of gauge invariant operators. In fact the gauge symmetry is not a symmetry of the dual theory in accordance with the usual understanding that a gauge symmetry is a manifestation of redundant degrees of freedoms that have to be eliminated. Taken to its logical conclusion, the principle of gauge invariance cannot be used as the basis to construct such a theory and perhaps more ideas are needed to understand the fundamental principles lying behind gauge theories.
For these reasons it is natural to study gauge invariant formulations of gauge theories. In fact it is known that gauge theories can be formulated entirely in terms of Wilson loops. In particular, in the large-N limit Wilson loops obey an equation that closes in the expectation value of single Wilson loops. This is known as the loop equation [3,4] or the Migdal-Makeenko equation. For finite N the equation is also valid but closes in the expectation value of disconnected (i.e. multitrace) Wilson loops. Since it is not clear how to renormalize the loop equation it can be better study perturbatively or in the lattice. Motivated by this, here we discuss different numerical and analytical methods to study the loop equation in the lattice.
Although, as we argued, this equation is of great importance there does not seem to be many studies on how to solve it. A notable exception is the very interesting work by Marchesini [5] where a formal solution and a numerical approach to solve the loop equation was proposed. We discuss this approach in detail but unfortunately it seems restricted to the strongly coupled regime, moreover, and the known numerical results are for the 2d case. The method we propose is based on the observation that certain matrices constructed of Wilson loop expectation values are positive definite. They can be thought of as reduced density matrices in the space of open loops after tracing over color indices. Imposing this extra condition allow us to select valid solutions of the loop equations. We call this approach a bootstrap approach since it uses general positivity properties of the theory to impose bounds on the solutions and also since imposing positive definiteness is an important part of the recently developed and highly successful conformal bootstrap program [6,7] 1 . The bounds we impose are for the expectation value of the energy. In two dimensions such bounds constrain the solution to be equal to the exact solution with high degree of accuracy. In four dimensions the bounds are less restrictive (due to larger computational complexity) but we can resort to a simple approximation, at small coupling we minimize the expectation value of the energy subject to the constraints, whereas at large coupling the entropy should be maximized. This gives results which are in good agreement with simulations and a reasonable approximation to the coupling where the transition occurs. The weak coupling region is well described by this method although, at the moment, this is not enough to understand the continuum limit.
For future work, it seems of great interest to apply this method to N = 4 SYM in order to make contact with the AdS/CFT correspondence. In this paper we take a few initial steps in this direction by briefly considering the bosonic sector of N = 4 SYM but leave a detailed study for future work.
It is interesting to note that the relevance of the extra positivity conditions at weak coupling was already observed [9] in the collective field method of Jevicki and Sakita [10]. Using the Kogut-Susskind approach a Lattice Hamiltonian in loop space was derived and numerically studied in [9]. See also the related work by Yaffe [11] using coherent states. This paper is organized as follows, in the next section we review the derivation of the loop equation and summarize the main ideas presented in this paper, following that we consider in great detail the two dimensional case since its solution is known exactly and can be used to test various approaches very easily. Afterwards, we apply those ideas to the four dimensional case and show that our numerical approach gives a good understanding of the gauge theory in the small coupling regime relevant for the continuum limit. Finally we describe a numerical simulation used to validate the results, discuss briefly the case of N = 4 SYM and conclude with a summary of the results and possible extensions and improvements.
2 Lattice gauge theory, a brief summary.
In this section we consider a four dimensional SU (N ) gauge theory in an infinite cubic lattice with Wilson action and briefly review known results for the large N limit. Then we discuss the derivation of the loop equation and introduce the notation we use in this paper. There is an extensive literature on the subject, our presentation here is just to summarize known results that are needed later in the paper and mostly follow the classic review [12] and the book [13] as regards to the loop equation.

Lattice action and known results
The system we consider is a cubic lattice where to each oriented link is associated a matrix U µ ∈ SU (N ). To the same link with opposite orientation we associate the matrix Uμ = U † µ . The action is the Wilson action where the sum is over all oriented plaquettes P . Here U P is the product of the four matrices associated with the plaquette and oriented means that we sum the trace of both possible orientation so that the action is real. The partition function is In four dimensions, for N ≥ 4 numerical results indicate that this theory has a first order phase transition as a function of λ [12]. In the large-N limit the transition occurs at λ c = 1.3904, as computed using the Twisted Eguchi-Kawai (TEK) model [14]. The nature of the transition is easily understood by considering the partition function in eq.(2.2) as defining a classical four dimensional statistical system with Hamiltonian and temperature T = λ. At small temperature λ → 0 we minimize the energy, i.e. the links U µ fluctuate around gauge trivial configurations. At large temperature λ → ∞ the entropy should be maximized and the links variables U µ explore all possible values with equal probability. Thus, the transition is a typical first order first transition between ordered and disordered states. The large coupling phase is confining and easily studied analytically in terms of a strong coupling expansion already proposed by Wilson [15]. On the other hand the continuum limit is obtained in the region λ → 0 which is more difficult to study. A Wilson loop expectation value is a real number associated with a closed path in the lattice and defined as where inside the trace we multiplied in cyclic order all the matrices associated to the given path C = {µ 1 µ 2 . . . µ L }.

Strong coupling phase
Analytically, the Wilson loop expectation values W C can be compute in a strong coupling expansion λ 1 by expanding the exponential of the action. The result has an interesting interpretation in terms of a sum over surfaces ending on the loop [16]. In 4 dimensions and in the large N-limit the expectation value for the plaquette is given by where we called the plaquette as Wilson loop one W 1 also denoted as u.
Wilson loop zero is the single point or null loop that obeys W 0 = 1. Such large orders in perturbation theory can be computed by using a character expansion of the exponential. Notice that W 1 = u is of particular importance since it determines the average Energy density (energy per lattice site) where V is the number of sites and the factor six is the number of plaquettes per site computed as d(d − 1)/2, in dimension d = 4.

Weak coupling
For small coupling the U µ matrices only have small fluctuations around the identity (up to gauge transformations). In that case one can write the theory in terms of a hermitian gauge field A µ , U µ = e iAµ and use Feynman diagrams to compute Wilson loops. The result for the plaquette in four dimensions and the large N limit is [17] Of course perturbation becomes unreliable for loops of large area since it does not capture the phenomenon of confinement that implies that such loops obey the area law.

Transition in mean field approximation
The transition from weak to the strong coupling regime is a first order transition that can be understood by a simple mean field approximation [12]. Within this approximation and in axial gauge, the free energy per site is given by [12] where v is a free parameter that is fixed by minimizing a, i.e. ∂a ∂v = 0. The function a(v) has a minimum at v = 0 and for λ < λ M 1.68 has a second minimum corresponding to the small coupling phase that has lower free energy for λ < λ c 1.48. The expectation value of the plaquette is given by as can be derived from the partition function. In the strong coupling regime we get the very crude approximation u = 0, in fact all Wilson loops vanish. In the small coupling regime we get a non-trivial function that can be expanded around λ = 0 as This agrees with eq.(2.7) but the λ 2 term (− 7 288 λ 2 ) is already incorrect. The mean field approximation can be improved [12] but we just wanted to emphasize that a simple approach captures the important physics.

Numerical simulations
All previous results are for infinite lattices and in the strict N → ∞ limit. For finite N and small lattices one can use numerical simulations and extrapolate the results to infinite N. The first results in this direction were by Creutz and Moriarty [18] and more recently by Meyer and Teper [19]. As part of this work we performed a numerical simulation for N = 10 in an 8 4 lattice which allowed us to check various ideas regarding the loop equation. The results for the expectation value of the plaquette are displayed in fig.1 where a very good agreement is seen with the perturbative and large coupling expansion in their regimes of validity. The position of the phase transition for SU (10) is seen to be around λ c (N = 10) 1.46 in agreement with the literature and close to the large-N value λ c 1.3904 [14].

Summary
Clearly there is already a very good and detailed understanding of this system, the strong and weak coupling regimes can be understood by series expansion and the transition using mean field. Numerical simulations validate the whole picture as summarized in fig.1. It should be emphasized however that the main physical interest lies in the continuum limit that appears in the λ → 0 region. In practice one has to show that large loops obey the area law in the weak coupling phase, a result that cannot be obtained by perturbation theory and can only be found by extrapolation of numerical simulations. In any case, our intention here is to study this system purely in terms of gauge invariant operators, namely with no reference to the variables U µ . For that reason we review now the derivation of the loop equation.

The loop equation
The loop equation is a direct consequence of the Schwinger-Dyson equation associated with a link of the lattice [3,4,13]. We reproduce its derivation here, first to introduce the notation, and second because for other theories the derivation will be done just by analogy to this one. Consider then a point x in the lattice and a given link µ = 0, 1, 2, 3, and perform the following change of variables The variation of the action is (2.14) iW ac where the terms in the last line come from the possibility that the path goes through the same link again either in the same direction (n + times), or opposite direction (n − times). These terms are call self-intersection terms but notice that self-intersection in this context means that the loop goes through the same link more than once (in either direction) and not merely through the same vertex. Since this identity is valid for any traceless hermitian we conclude (2.17) Contracting both sides with δ ac δ bd we get and used the large N factorization property The different terms in the result have a very simple graphical interpretation as seen in figs. (2,3). The link x, µ appears in the action in several plaquettes. Each of these plaquettes is connected to the loop as in the figure, when the orientations are opposite, we include a minus sign. Also the action comes with a coefficient N 2λ . Summing over all links we get the loop equation that we schematically write as where L denotes the length of the loop, S * W denotes all possible intersections between the loops appearing in the Wilson loop (at a fixed position) and those appearing in the action. The last term is a sum over all self-intersections with a sign σ i depending on the orientation of the intersection and C 1 and C 2 denote the two loops in which the original loop splits when reconnecting at that selfintersection (see fig.3). In this form the loop equation is valid for any action given by a sum of Wilson loops. Therefore, we can give a linear combination of loops as the action and the reconnection procedure determines completely the loop equation without any reference to matrices, gauge invariance etc. For mathematical manipulations it is convenient to enumerate the loops in a list, where we eliminate redundancy due to rotations, translations, cyclic permutations and opposite orientations. The first few elements of the list are in fig.(4). Then we can write the loop equation in the form where W i denotes the expectation value of loop i, K i→j is a matrix indicating that loop i converts into loop j by the reconnection procedure of the action with a weight depending on how many different ways we can get j and divided by the length of the loop L. The tensor C i→jk is the self-intersection term and indicates that Wilson loop i splits into jk with an appropriate coefficient. For example, for the plaquette we get

Extra equations
From the derivation of the previous subsection it is clear that we can get more equations than just the loop equation. First we can obtain individual Schwinger-Dyson equations for each link. Given two links that do not belong to a self-intersection, the difference between their respective equations is linear and independent of λ. Other linear lambda-independent equations can be obtained from the equations associated to links that touch the loop but do not belong to it (see fig.5). All these equations are linear and λ-independent and we denote them as constraints since do not have information on the coupling: An example is: They should also be imposed since they restrict the possible values of the Wilson loops.

Two dimensional lattice
The two dimensional system is a well known system that can be solved exactly even in the large-N limit [20]. In this paper we use it to test numerical methods that can then be extended to the more challenging case of four 0 Figure  dimensions. In axial gauge, U 0 = I, the two dimensional case reduces to the single plaquette: The Wilson loops can be labeled by an integer The large N limit was studied by Gross and Witten as well as Wadia [20] using the saddle point of an effective action for ρ(θ), the density of eigenvalues e iθ of U in the interval θ ∈ [−π, π]. Later Friedan [22] obtained the same result using the loop equation. The result for the plaquette (and therefore the energy) is At λ = 1 there is a jump in the second derivative of W 1 and therefore the transition is third order. Since the solution is exact we can use this system to test different methods to solve the loop equations as we do in the rest of this section. In terms of the eigenvalue density ρ(θ) (normalized to π −π ρ(θ)dθ = 2π) where we used that the W n are real. Thus, the eigenvalue density is a generating function for the Wilson loops.

Positivity constraints, density matrix and Wilson loop entropy
The fact that the matrices U are unitary imply certain constraints that are fundamental to understand the physics and to implement the numerical methods that we describe below. We start with the observation that for any matrix A with components a ij : whereρ 0 indicates thatρ is positive semi-definite. Since W 0 = 1 then Trρ = 1 and thereforeρ has properties of a density matrix. Indeed, we can write its definition asρ where U ( ) = U . Namely, if we take the collection of all powers of the matrix U up to power L, the matrixρ traces over the color indices, the entropy S W L = −Trρ (L) lnρ (L) measures the information loss due to such tracing. Mathematically, the matrixρ (L) is a Toeplitz matrix defined as in the eq.(3.7) or, equivalently, A useful comment is that the Toeplitz matrix is associated with the Fourier coefficients of the eigenvalue density ρ(θ) of U . In such case one can use Szegö theorems [23] to compute limits of functions of the eigenvalues ofρ (L) by using integrals of the eigenvalue density: j=1...L are the eigenvalues of ρ (L) and ρ(θ) is the eigenvalue density 3.4.
From these constraints one can derive some simple results that are completely independent of the action that we choose: • For example taking the principal minor implies that all loops satisfy |W n | ≤ 1 as we already know. Taking another principal minor we obtain and therefore, if W 1 = u = 1 then all loops are equal W n = 1 independently of the action! If u < 1 we have an interesting inequality that bounds the rate of change in the Wilson loop expectation value We can then look for actions that saturate these bounds, an interesting topic that we leave for future work.
• Notice, from the previous item, that the bound forρ (L=2) is |u| ≤ 1 and when we saturate it (u = 1) we are able to compute all Wilson loops obtaining W n = 1. This is generic. Suppose we saturate the bound for ρ (L) . That happens when we develop a zero eigenvalue ofρ (L) , namely there is a particular set of coefficientsĉ i such that But this is the mean value of non-negative quantities and therefore can vanish only if it vanishes for all configurations: for a specific set of coefficients {ĉ i=0...L }. Notice this is a matrix equality that means we are considering only configuration that satisfy this specific constraint. Now compute valid for any j ∈ Z. Take the largest r such thatĉ r = 0 (normally r = L) then which is a recursion relation that allows us to compute all Wilson loops assuming that we already know the ones up to W L . Again it should be interesting to find models that respect these equations. In the case of the Wilson action equations such as this are not exact. Nevertheless we expect a similar equation to be valid in the limit L → ∞.
• The space of positive definite matrices is a multi-faceted convex cone. One can envision phase transitions when the minimum of the action jumps from the interior to the boundary of the cone, or between faces in the boundary. As we see later, for the case in hand the transition is of the former type. For finite L, the strong coupling solutions lies in the interior and the weak coupling ones at the boundary.
• The matrixρ is positive definite for any value of N and λ. For example for N = 1 the Wilson loops are given simply by Bessel functions implying that the matrix T[I 0 (1/λ), . . . , I L (1/λ)] is positive semi-definite for any L and λ.
We emphasize the previous points since they translate also to higher dimensions as we discuss later.

Effective action, numerical solution
In [20], the following effective action for the eigenvalue density was constructed we get, up to an additive constant, a simple effective action for the Wilson loops Minimizing this action with respect to the W n trivially gives W 1 = 1 2λ , W n≥2 = 0, namely the strong coupling solution. However, for λ < 1 the corresponding matrices ρ L are not all positive definite. Indeed, their eigenvalues are µ . . L + 1 that are not all positive for λ < 1. Therefore the correct problem to solve is for some fixed L. This problem has the form known as Quadratic Programming and can be converted into a problem of Semi-Definite Programming (SDP) (see [24] and the appendix) and solved by standard packages. Using the matlab cvx package [25] it is very easy to show numerically that for L = 10 the results for W 1 = u(λ) agree perfectly with the exact answer, and even for as low as L = 6 they are reasonably correct. Increasing further L one can get more precise values for W 1 and also compute the other loops W n=1...L , in good agreement with the exact answer. This approach provides an excellent numerical method purely in terms of the

The Loop equation and exact solution
Consider the loop equation for a given loop W n>0 of length L = 4n. Every link gives rise to the same result, when connecting with a plaquette we get W n+1 − W n−1 . Each link self-intersects with n − 1 other links giving W p W n−p for p = 1 . . . n − 1. After dividing by the length, the loop equations becomes simply [22,13]: with the condition W 0 = 1. It is clear that, if we give a value to W 1 = u, then all the other Wilson loops are fixed recursively. The equation is very powerful but unfortunately it still leaves an infinite number of solutions, one for each value of u. The first observation is that |W n | ≤ 1 since |TrU n | ≤ N because U is unitary. A little numerical experimentation shows that, for λ > 1 the recursion leads to divergent values of W n as n grows, except if u = 1 2λ . If λ < 1, the recurrence diverges if u < 1 − λ 4 but is finite for u ≥ 1 − λ 4 so more constraints are needed at small coupling. Formally, following [22], we can define a generating function The condition that |W n | ≤ 1 implies that ϕ(z) is analytic inside the unit circle. This fixes u = 1 2λ at strong coupling but allows any 1 − 1 2 λ ≤ u ≤ 1 at weak coupling. The extra condition that we need at weak coupling is that the eigenvalue density is non-negative which is enough to determine ϕ as shown in [22]. The reason we repeat it here is that we wanted to emphasize the main message: • There is an infinite number of solutions to the loop equation.
• The constraint |W n | ≤ 1 greatly reduces the set of solutions, especially at strong coupling.
• An extra condition such as eq.(3.28) is needed in order to find a unique solution.
These properties can be translated into a simple numerical method which can be extended to four dimensions. Before describing it, we discuss a previous method due to Marchesini [5] that produces good results at strong coupling but not at small coupling emphasizing the difficulties encountered in that region.

Numerical solution: Marchesini's approach
A simple numerical method to solve the loop equation was proposed and tested in two dimensions in [5]. Since the number of Wilson loops is infinite, any numerical approach has to chop the set of loops. Let us assume that we consider loops up to length 4L, namely W L is the last one. Given The roots of the polynomial give the possible values of u. This polynomial has always a root u = 1 2λ that corresponds to the strong coupling solution. For small λ other solutions appear. For example, for L = 49 we plot the roots in fig.6. It is clear that the small coupling solution will appear in the limit of large L as an envelope of the lowest roots. However, the figure does not allow us to expect a nice convergence. In [5] an iterative method is proposed that converges to the lowest root and therefore it should give the correct value of u. However, in the same reference it is pointed out that the method requires various cut-offs and extrapolations in the weak coupling region. This suggest that its four dimensional formulation might be hard to deal with. Now we turn to a simple method that we propose in which a different polynomial whose roots provide upper and lower bounds to u that quickly converges to the actual value (already L = 8 matches the exact solution).

Numerical approach: Bootstrap-like approach
The main point of the numerical method is to implement as many constraints from unitarity as possible. After fixing a maximum L, instead of setting W L+1 = 0 we allow it to vary under the constraint thatρ (L) 0 and thus find analytical bounds for the expectation value of W 1 = u.
Before doing that, however, we implement an even simpler idea that is to impose just the constraints |W n=1...L | ≤ 1. Starting from L = 2 and incrementing L, for fixed L, as we vary u the first loop to violate the constraint is the largest one, W L . Therefore this is equivalent to set W L (u) = ±1. From the roots of this polynomial in the region allowed by the previous step L − 1, we choose the smallest one. This already gives better results at small coupling as can be seen in fig.7. The curves are lower bounds that are continuous and converge to the exact solution as L is increased (we reach L = 36).
Imposing all the constraints contained inρ (L) 0 is even better. In the allowed region,ρ 0, all eigenvalues ofρ (L) are positive. As we vary u, we reach the boundary of the allowed region when an eigenvalue vanishes, namely detρ (L) = 0. This determinant is a polynomial in u and therefore the roots of such polynomial determine the analytical bounds of the allowed region, for given L. Again we increase L by one in each step contracting the allowed region every time. The results are displayed in fig.8 where we can see that already for lower values of L the approximation is very good. For reference we give The two roots of this polynomial contained in the interval [0, 1] correspond to the blue curves in fig.8.

Small coupling expansion
At small coupling we can use the approximation   vanishing term gives: which then implies w 1 = 1 2 as we know from the exact result. Higher orders should vanish as can be obtained by considering larger values of L. It is absolutely clear then that a systematic expansion in λ for small coupling can only be achieved by using theρ matrices. The loop equation alone is not enough to fix this expansion.

Approximation
Although the results are excellent already for small values of L we can consider a relatively low value, e.g. L = 4 and wonder if it is possible find an approximate value of u between the maximum and the minimum. This would be an approximation as opposed to the bounds that are analytical bounds. The low coupling phase is a low temperature and therefore should minimize the energy, equivalently maximize u. So, at small coupling we choose u = u max which indeed gives a very good answer, see fig.8. The large coupling or large temperature phase should have large entropy. In this case there is a simple entropy we can consider, namely the entropy associated to the matrixρ (L) . Indeed, in the limit λ → 0, all loops are given by W n = 1, ρ (L) M has one eigenvalue equal to one and all the others vanish. Namely it describes a pure state. Up to gauge transformations, the matrix U is the identity and we do not lose any information if we trace its powers. On the other hand when λ → ∞ all loops vanish except W 0 = 1. The matrixρ is proportional to the identity and the entropy is a maximum. Namely, if we take traces of powers of U we lose a maximum amount of information for these configurations.
Therefore, the simple proposal is to maximize u for small coupling and maximize S W L at large coupling. The results agrees with the exact solution better than the bounds. An even better result is obtained by defining an effective free energy where c is an adjustable constant of order one. We set c = 1 2 because it seems to adjust the exact answer well (see fig.9) but we do not have a way to fix this constant from first principles. Of course the correct effective action is the one we gave in section 3.2 but here we wanted to find a simple effective action that could be used also in higher dimensions.

Summary
To summarize, what we learned from the simple two dimensional case is: There is an infinite number of solutions to the loop equation but they are restricted by imposing positivity ofρ (L) . Such condition also allows the derivation of bounds independently of the action. Once the loop equation is imposed we can derive a weak coupling expansion, strict bounds on the energy and a simple approximation when considering short loops.

Four dimensional lattice
In four dimensions the numerical methods are similar as in two dimensions, the main difficulty being that the number of Wilson loops grows exponentially with the length. To handle that, we developed a computer program that listed all loops up to translations, rotations, reflections and cyclic permutations of the links up to length L = 18 although most calculations described below were done using loops up to length L = 14. It also computes the corresponding loop equations, strong coupling expansion, and a set ofρ matrices that have to be positive definite as explained below. Finally it provides output that can be further manipulated by computer algebra programs or standard packages such as cvx [25] or sdpa [26]. Let us now briefly describe different methods that can be used to solve the loop equation in different regimes and their usefulness.

Strong coupling expansion for plaquette expectation value
The strong coupling expansion for the Wilson loop can be done straightforwardly using the loop equation [5]. Indeed writing the loop equation as we obtain a simple solution as a series expansion We obtain the expansion for the plaquette as Up to the computed order, the result agrees with eq.(2.5) thus providing a way to validate our computer code. Higher order terms require going to larger loops or using other methods such as character expansion [12].

Iterative strong coupling numerical solution
Instead of doing an analytical expansion we can do a simple numerical iteration of eq.(4.1).
As seen in fig.10, the results match very well the strong coupling expansion but diverge for λ 1.5. The reason is that, for the iterations to converge, the eigenvalues of the operator − 1 2 K have to have modulus less than λ. These eigenvalues are plotted in fig.11 where one can see that indeed the largest eigenvalue is λ max 1.5.

Strong coupling solution, analytic continuation to weak coupling
We can find a different iteration method to solve eq.(2.23) by doing (4.7) Clearly the iteration is ill-defined if λ is an eigenvalue of − 1 2 K. The matrix K is not symmetric but numerically we can still diagonalize it after truncating it by setting to zero Wilson loops larger than a certain length. For example, keeping loops up to length L = 10, we observe that the eigenvalues of − 1 2 K approximately cover an interval λ ∈ [−1.5, 1.5] on the real axis plus some sporadic eigenvalues in the complex plane that are likely the result of the truncation. We expect that in the limit of infinite length there is a cut on the real axis. By taking complex values of λ one can find a simple analytic continuation of the string coupling solution to small values of |λ| away from the real axis. However, one can see that such analytic continuation is not the small coupling solution as one can actually expect on general grounds since the transition is first order. This method also shows that the divergences on the real axis are due to the fact that we truncated the loops by putting the larger ones to zero. If that were not the case we could adjust them so that the right hand side of eq.(4.7) does not contain the problematic eigenvectors of − 1 2 K thus avoiding the divergences. The question arises of how should one choose the higher loops. This takes us to the next method, namely we choose them so that that matricesρ (L) are positive definite.

Positivity constraints
We have to construct the analogue of the matrixρ (L) in two dimensions. The idea is very simple, take two points  In this way we construct a set of matrices that have to be positive definite. Before going into the numerical procedure let us notice a few facts regarding this matrices that are independent of the action that we use to average the Wilson loops.
• If we take x 1 = x 2 and two loops that start and end at x 1 and share the first link, see fig.12, we can construct the matrix • Another example is a long loop a made out of two plaquettes connected by a long path ( fig.12). Cutting this loop in two as in the figure and using the same procedure we obtain where u is the expectation value of the plaquette. This means that, if the expectation value of the plaquette is close to one, there are arbitrary long loops that are also close to one.
• Finally notice that, as in 2D, if the matrixρ has a zero eigenvalue, namely it is in the boundary of the allowed region, then for a given set of coefficientsĉ in eq.(4.8) the matrix A vanishes and this in turn implies an infinite set of linear equations for loops that contain the given open paths. Namely, taking an arbitrary path C 0 from x j to x i we get Tr U 0 U ( ) = 0 . (4.12) As a simple example, let us show that if the plaquette is one (u = 1) then all loops are equal to one. Indeed, take two paths that build a plaquette as in fig.13, we get the matrix where u is the plaquette, and the third one implies that if u = 1 then for all configurations the matrices U 1,2 associated with the two paths are equal U 1 = U 2 and then all loops are equal to one.
as expected. However, if the bound is saturated, i.e. u = 1 then the matrix associated with the zero eigenvalue should vanish. Namely A = U 1 − U 2 = 0. Since we can use this in any loop it simply means that, for any loop, performing a move such as the one in the same figure does not change its expectation value. Since any loop can be reduced to a point by such moves, then all loops are equal to one.

Bootstrap-like method
The numerical method should now be clear. We list all loop equations and extra equations we described previously and that can be handled by the computer resources available. Then we construct all matricesρ that we can handle and look for solutions of the loop equation that satisfy the constraint ρ 0. This gives an upper and lower bound on u. Let us first do this analytically for the simplest example, the loop equation for the plaquette in dimension d ≥ 3. In the notation of fig.4, the equation is (4.14) If we want to maximize u, we set W 2 = W 3 = 1, their maximum values and W 20 = W 17 = W 21 = 2u 2 − 1 their minimum values according to the previous subsection, eqs.(4.10) or (4.11). We then get an equation for the maximum One of the roots of this equation gives the upper bound, namely a bound valid for the Wilson action in a cubic lattice of any dimension d ≥ 3. For λ → ∞ we get u max 2(d−1) λ and for λ → 0, u max 1 − λ 4(d−1) . The bound has the right behavior but the coefficients are not right, as was somewhat expected from such crude bound. This calculation is simply an illustration that there is indeed a bound that follows from positivity ofρ and the loop equation. Similarly we can get a lower bound by choosing W 2 = W 3 = 2u 2 − 1, and W 20 = W 17 = W 21 = 1: (4.17) To go further we can use a numerical procedure that allows us to handleρ matrices of size of order 10 3 × 10 3 and thousands of equations. The main obstacle is that the known numerical packages (see appendix) only deal with linear equations. The loop equation is non-linear but we can go around it by fixing a set of loops such that the equations become linear in the rest of the variables and then exploring the available space of such loops. In this paper we only consider the case where we fix the plaquette u, and therefore we have to explore only a one-dimensional space thus simplifying the calculation. In practice we consider loops up to a given maximum length L and the actual procedure depends on L. Let us now consider each case.

Linear case, L max = 8, 10
Since we consider Wilson loops of maximum length given by L max = 8, 10 we can only impose the loop equation associated with loops up to length L = 6. Those loops do not have self intersections and therefore the loop equations are linear. In this case we can solve the problem using semi-definite programming in a direct way. For example using matlab and cvx [25] (see appendix) we find the bound depicted in fig.14. We already see that the results are reasonable for the maximum value of u at small coupling. The minimum value of u for L 8,10 turns out to be u min = 0 which is a correct but rather poor bound for the actual value of u. So, we consider now Wilson loops of larger length.

Non-linear case in one variable, L max = 12, 14
In this case we impose the loop equation up to length L = 10. Some of those loops self-intersect and we cannot use semi-definite programming directly. However, the self-intersection splits the loops into two whose total length is less or equal than L = 10 and therefore one of them at least has to be a plaquette. For that reason we propose a different semi-definite programming problem. We fix the value of u and define a new matrixρ aŝ Now we maximize the value of t by allowing the loops other than the plaquette to take arbitrary values compatible with the loop equation. When the procedure finalizes, the value of t is equal to the lowest eigenvalue ofρ and it is the largest lowest eigenvalue that can be found for that fixed value of u. Thus, if t max < 0 it is simply not possible to choose the other loops such that the loop equation is satisfied and the matrixρ 0. Therefore this value of u is not allowed. In that way we can sweep the allowed values of u and determine the bounds on u and therefore on the energy. The results are depicted in fig.14 where the matrixρ was truncated to an 800×800 size. The bounds are not close to each other as they were in two dimensions. At small coupling we minimize the energy and therefore choose the maximum value of u. Following the ideas discussed in section 3.7 for the two dimensional case, for large coupling we should maximize the entropy of the matrixρ. This is not an SDP problem and therefore we do a further approximation. Consider the value of u where t max (u) has a maximum. One can associate such point with a large value of entropy since increasing the minimum eigenvalue, with the trace being fixed, tends to make all eigenvalues similar. We are going to take such value of u as the best guess of the one that maximizes the entropy. This is depicted in fig.14 where we see that it is quite a good approximation at strong coupling.

Small coupling expansion
As in section 3.6, we can do a small coupling expansion by linearly expanding all loops as and finding bounds on w n . Using sdpa for loops up to length L = 12 we find for the plaquette w 1 ≥ 0.238 and for L = 14 w 1 ≥ 0.2495 in agreement with the value w 1 = 1 4 . Unfortunately for this length we do not find a minimum value of u and we cannot fix w 1 = 1 4 as we were able to do in two dimensions.

Lattice simulation
Our Monte Carlo simulation used the Metropolis Algorithm to produce an ensemble of uncorrelated configurations using the partition function in eqs.(2.1), (2.2): (4.20) In order to ensure that the results were properly thermalized we allowed the simulation to run for 20000 updates before saving configurations. Each update consisted of 10 "hits" on each link. A hit is one attempt to move the link through the phase space. After thermalization, one configuration every hundred was saved until a total of 300 configurations were collected. Binning showed that the correlation time of the system was around 400 updates between saves. Error was calculated via the bootstrap method and we found that for Wilson Loops up to length 10 the error in the expectation value was around 10 −5 for lattices of size 8 4 . The simulation was programmed using CUDA on GeForce GTX 980s.  fig.15 the results for loops W 2 , W 3 , W 4 (see fig.4) are displayed. The agreement is reasonable at small coupling λ 0.5 if we used the solution that minimizes the Energy (maximizes u), in agreement with our previous discussion. For larger values of the coupling but below the transition the agreement is qualitative in the sense that the larger loops have smaller values than the plaquette. In the strong coupling phase, we choose the solution that maximizes t but the previously found agreement ( fig.14) for the plaquette does not extend to the other loops. Clearly, the well-known strong coupling expansion is still the best method in this region.

N = 4 SYM
The case of N = 4 SYM is particularly important since it has a dual description as a string theory [1]. In the context of its relation to string theory and more precisely with the AdS/CFT correspondence, the loop equation has been formulated and studied in [28]. Studying the theory in the lattice should allow a different method of computation and possible strong coupling calculations based on the gauge theory side of the correspondence. In the rest of this section we briefly describe how the ideas we developed for pure YM could be implemented but, for concrete calculations, at the moment we restrict ourselves to the bosonic sector leaving the study of the full theory for future work.
A lattice theory that has the correct continuum theory without the need for fine-tuning is described in [29]. Here we use that formulation but follow the notation found in the paper [30]. For brevity we do not explain details and refer the reader to that work for explanation of the notation and properties of the theory. The main property is that such formulation preserves one twisted [31] scalar supercharge and possesses BPS Wilson loops although more restricted than the continuum theory. For our purposes another important property is that the action can be formulated entirely in terms of generalized Wilson loops, namely using loops with fermionic links and/or sites. In this way, the form of the loop equation given in (2.22) is valid using the appropriate supersymmetric action [30]. However, we have not worked out the correct generalization of theρ matrices to the fermionic sector and therefore here we restrict ourselves to the bosonic sector described by the simpler action This action, thought as a linear combination of Wilson loops is depicted in fig.16. In order to obtain the loop equations we must vary individual links. However, the scalar and the gauge fields are twisted together which means that the matrices U a (x) ∈ GL(N, C) and therefore the links and their daggers must be treated independently: This implies that in eq.(2.22), the intersection of the action with the loop is non-zero only for links oriented in the same direction, the same is true for a loop self-intersection. Since the gauge group is U (N )

Monte Carlo Simulation
We used the parallel code developed by Schaich and DeGrand [32], based on the previous work by Catterall and Joseph [33] to simulate the latticized Euclidean theory. The code allows for a simple way to reduce to the bosonic sector, i.e. to the six scalars and four gauge fields that in this formulation live on the links. The lattice is taken to be an A * 4 lattice which contains the permutation group S 5 , the largest finite subgroup of the four dimensional rotation symmetry. In order to preserve a subgroup of the SUSY Algebra, the Euclidean-Lorentz and R symmetry groups are twisted into SO(4) E ⊗ SO(6) R → SO(4) ⊗ U (1). With this twisting the lattice theory is invariant under one supercharge out of the full sixteen. The continuum limit should restore the full sixteen supercharges without fine-tuning.
To test the loop equations we considered the loop equation associated with the plaquette and found that it is satisfied up to four digits which is within the numerical accuracy of the simulation. The test was done for U(2) and U(3) gauge groups and for λ lat = 0.

Conclusions
The loop equation has been traditionally seen as a promising way to describe gauge theories in terms of gauge invariant quantities. In this paper we agree with this perspective but also point out that such equation has infinite solutions that have to be constrained by the condition that certain matrices are positive definite. At strong coupling such extra conditions are not necessary and therefore seem to have been largely ignored. On the other hand, in the physically relevant region of small coupling such conditions are crucial to obtain the correct solution. In fact they also give many constraints and properties that are actually independent of the action. In two dimensions this leads to a simple numerical procedure that reproduces the exact solution for any coupling. This method can potentially be used for other two dimensional actions where the exact solution is not known but, more importantly, it can be extended to higher dimensions. The simple idea is that given two points in space and a set of open lines = 1 . . . L between them, we can define an L × L matrix of closed loops where the entry is the expectation value of the closed loop obtained by going along path and returning along path . Such matrix has to be positive definite and can be considered as a reduced density matrix due to tracing over the color indices. Its entropy measures how much information is lost by taking the traces and gives a qualitative idea of the entropy of the system. The reason is that the entropy of the matrix ρ vanishes in the ground state when all links are equal to the identity (up to gauge transformations) and is maximal at large coupling when the links fluctuate randomly. Numerically we implement a procedure to solve the loop equation under the condition that the constraints are satisfied. It reproduces known results from Monte Carlo simulations but cannot be considered an improvement. It is possible that more computational resources could lead to better results especially in the small coupling region relevant to the continuum limit. Also, it might be possible to improve the choice of the constraints, namely the matrixρ, and/or use improved actions that approach the continuum limit faster. It is also interesting to extend this method to other theories whose action can be formulated entirely in terms of closed loops. One such theory is N = 4 SYM, of particular importance since it plays a central role in the AdS/CFT correspondence. Initial steps in this direction suggest that the method can be applied but requires a better understanding of the fermionic loops. Another system to consider is three dimensional Yang-Mills where there are other approaches [35,34].

Acknowledgments
This work was supported in part by the DOE through grant DE-SC0007884. We are very grateful for numerous discussions with P. Vieira, S. Catterall, D. Schaich on the matters of this work and/or lattice gauge theory in general. Also S.Catterall and D. Schaich graciously provided their lattice code allowing us to check the loop equations in the bosonic sector of N = 4 SYM. We are also grateful to D. Minic and A. Jevicki for useful comments on a previous version of this work. P.D.A. would like to thank the Wigner GPU Laboratory at the Wigner Research Center for Physics (Budapest, Hungary) for providing GPUs computer resources. He would also like to thank G. G. Barnaföldi, M. F. Nagy-Egri, D. Berényi, and Z. Bajnok for helpful discussions. M.K. wants to thank the hospitality of the Perimeter Institute, Waterloo, CA and the SAIFR Institute (São Paulo, Brazil), while part of this work was being done.

Appendix: Semi-Definite programming (SDP)
Semi-definite programming [24] is a type of optimization problem that has been the focus of a lot of attention recently in relation to problems in finance, engineering, and more recently has proved invaluable in the bootstrap program of conformal field theories (see e.g. [7]).
It can be stated very simply as: Given m real numbers c i=1...m ∈ R and m real, symmetric n × n matrices F i=0...m ∈ R n×n , find x i=1...m ∈ R that minimize m i=1 c i x i under the constraint that X = m i=1 x i F i − F 0 is positive semi-definite (X 0). The main observation in this field is that the space of semi-definite matrices is a convex cone in the space of all symmetric n×n matrices. Indeed, given two positive semi-definite matrices X 1,2 0, namely y t X 1,2 y ≥ 0, ∀ y ∈ R n , then it is clear that y t (α 1 X 1 + α 2 X 2 )y ≥ 0 for any real α 1,2 ≥ 0. Thus, α 1 X 1 + α 2 X 2 0, ∀ α 1,2 ≥ 0 showing that positive semi-definite matrices form a convex cone. The condition that X belongs to the linear subspace generated by the F i=1...m shifted by F 0 defines an intersection between this hyperplane and the semi-definite cone. This is a convex region over which we minimize a linear function. Therefore, the minimum is unique and should be located at the boundary of the domain, namely when the matrix X has at least one zero eigenvalue. The problem is then very similar to the more tradi-tional problem of linear programming where one minimizes a linear function m i=1 c i x i over the convex cone y =1...n ≥ 0 where the y = m i=1 a i x i for some given coefficients α i .
There are many other problems that can be reduced to an SDP problem. For example, in section 3.2 we need to solve this is problem of quadratic programming that can be reduced to an SDP problem by defining a new variable t and imposing or equivalently  Once the problem has been casted as an SDP problem, there are several available packages that can be used to solve it. We found that for rapid development of small problems the matlab package cvx [25] was very convenient and, for larger problems sdpa [26] was a good choice.