Testing volume independence of SU(N) pure gauge theories at large N

In this paper we present our results concerning the dependence of Wilson loop expectation values on the size of the lattice and the rank of the SU(N) gauge group. This allows to test the claims about volume independence in the large N limit, and the crucial dependence on boundary conditions. Our highly precise results provide strong support for the validity of the twisted reduction mechanism and the TEK model, provided the fluxes are chosen within the appropriate domain.


Introduction
Gauge theories lie at the basis of our understanding of Particle Physics. They are also beautiful mathematical models which are full of rich phenomena, both perturbative and non-perturbative. The pure gluon version of the theory has no free parameters and its dynamics underlies many of the fascinating properties of the interaction of quarks, such as Confinement and chiral symmetry breaking. Despite having no scale in its formulation, the theory is not conformal invariant because a scale Λ is generated by the mechanism of dimensional transmutation [1]. All physical properties of the theory can then be expressed, depending on its dimensions, in the appropriate units given by powers of Λ. Observable quantities which involve large energy scales (with respect to Λ) can be calculated by perturbative methods, thanks to the property of asymptotic freedom [2]- [3]. On the contrary, if some low energy phenomena is involved, perturbation theory fails dramatically and one has to use non-perturbative techniques. Fortunately a very powerful formulation was introduced by K. Wilson [4], allowing the development of several non-perturbative methods of computation. Lattice gauge theories define the quantum field theory as the limit of an statistical mechanical system at critically. One can then employ numerical methods which were originally devised to handle statistical mechanical systems. Unfortunately, for the study of the critical behaviour only numerical Monte Carlo methods seem useful. Nonetheless, progress in computer technology and in numerical algorithms have made possible the calculation of several observables of Yang-Mills and the more demanding theories with dynamical fermions (for recent results see Ref. [5]).
In the 70's G. 't Hooft [6] suggested the use of a new expansion parameter: 1/N , the inverse of the rank of the group. Asymptotic freedom, Confinement, dimensional transmutation and a massive spectrum are still present at the lowest order approximation: the large N limit. Furthermore, it is known that several aspects of the theory are simpler in that limit: planar graphs, stable spectrum, factorisation [7], etc. It seems clear that a more profound understanding of the large N limit would precede that of finite N. The large N limit also plays a crucial simplificatory role in the new ideas and methods that have emerged under the heading of AdS/CFT [8].
One of the obvious questions is then whether the large N limit also introduces a simplification in combination with the lattice gauge theory formulation. At first, it seems just the opposite. Numerical methods are devised to work with finite N values. Results at the large N limit then follow from extrapolating those obtained at finite N . But, the numerical effort grows with the number of degrees of freedom (hence, at least, as N 2 ). Nevertheless, following this procedure might be worth the effort, in order to produce predictions that can be matched with those of other approaches. A large literature has been generated following this line which can be consulted in Ref. [9].
There is, however, one potentially simplificatory phenomenon which was discovered when combining large N with the lattice approach: reduction or volume independence. The idea emerged from a crucial observation made by Eguchi and Kawai(EK) [10] when looking at the loop equations [11] obeyed by Wilson loops on the lattice. Assuming invariance under Z 4 (N ) centre symmetry and factorisation, they concluded that these equations are independent of the lattice volume.
Thus, if the phenomenon is correct, one can reduce the study of the large N limit to a simple matrix model with only d matrices (d is the space-time dimension). This is, no doubt, a very powerful simplification, at least at the conceptual level. The first question is then whether the idea is indeed correct. However, there are other secondary questions which could be crucial for applications. In particular, one might inquire about the size and nature of the corrections. What are the leading corrections that one has at large but finite N? These are the central issues which our paper tries to address.
The history of the subject is long. The initial stages took place soon after the EK paper. It was shown that the EK model, obtained by reducing to a one-point lattice with periodic boundary conditions, violates the centre-symmetry condition for reduction at weak coupling [12]. Indeed, the one point lattice had been studied before the EK paper in Ref. [13], and the weak coupling path integral was shown to be dominated by symmetry breaking configurations. Fortunately, the authors of Ref. [12] proposed a modification under the name Quenched Eguchi-Kawai model (QEK) devised to restore the symmetry and, hence, the validity of reduction. On the other hand, the studies of Ref. [13] suggested that boundary conditions have a very strong influence on the weak coupling behaviour at finite volume. Thus, the present authors [14]- [15] presented a new version of reduction based upon twisted boundary conditions and aimed at respecting a sufficiently large subgroup of the centre symmetry group. The new model, called Twisted Eguchi-Kawai model (TEK), follows by imposing twisted boundary conditions and subsequently reducing the model to one point. At large volumes the boundary conditions do not influence the dynamics of the model, but at small volumes they are crucial.
Both the QEK and the TEK model passed all tests of absence of centre symmetry breaking done at the time of their proposal, and looked as valid implementations of the reduction idea. Unfortunately, attempts to solve these models explicitly failed. With an increasing interest in large N gauge theories the situation was re-examined in recent times, and signals of symmetry breaking were observed in both cases. By construction the QEK model is invariant under a certain subgroup of the Z 4 (N ) group. However, there are possible subgroups associated with the action on several directions simultaneously that can be violated and were observed to do so in Ref. [16]. In the case of the TEK model, although the remnant symmetry cannot be broken at infinite weak coupling, it can still be broken at intermediate values. The first observations of breaking were obtained by one of the present authors and collaborator [17]. Later on, several other authors [18]- [19] confirmed the breaking and showed how this problem affects crucially the possibility of reduction in the continuum limit.
During these years there have been other attempts to recover the validity of the reduction idea and its simplificatory character. Some authors argued that other large N gauge theories might be free from the symmetry breaking problems that invalidated the reduction of the EK model. In particular, theories with several flavours of fermions in the adjoint representation [20] are good candidates. This has triggered several efforts to obtain information about these theoretically interesting theories, by benefiting from the computational advantage of their reduced versions [21][22][23][24][25][26][27][28]. Although this is certainly very interesting, here we will not consider these theories any further since our interest in this paper is centred upon pure gauge theories.
Another approach was proposed by Narayanan and Neuberger [29], lately referred as partial reduction. The idea is to combine space and group degrees of freedom in an efficient way. Their study shows that at any value of the coupling there is a range of values of L (L being the linear size of the lattice box) for which the symmetry is respected at large N . More precisely, there is a minimum box size L c (b), which grows as the coupling gets weaker. Infinite volume large N results can be approached by tuning L and N appropriately. Although this idea was mostly developed with periodic boundary conditions, it is also possible to combine it with twisted boundary conditions [30]. This is a possibility that we will analyse further in the present paper.
Reasons for the failure of the TEK model were given in Ref. [17]- [18]. Starting from these studies, the present authors proposed a simple modification of the original procedure which could restore the symmetry and recover the validity of reduction [31]. The idea is to tune the fluxes appropriately when taking the large N limit. Indeed, the twisted boundary conditions are not unique; they depend upon the choice of discrete fluxes through each plane. In Ref. [15] we determined the range of fluxes that guarantee the preservation of the symmetry at very weak coupling. However, it is now clear that the actual choice of fluxes within this range has an influence at intermediate values of the coupling. In Ref. [31] we argued that, with a suitable choice, the problems can be avoided at all values of the coupling.
In the last few years we have set up a program to study this problem. As explained earlier, the question is two-fold: Is reduction correct? Is it useful? The second aspect might be crucial for its ultimate interest as a valid simplification. Trying to answer the second question directly, we set ourselves the goal of computing an observable of the large N continuum infinite volume theory using the TEK model. For simplicity we focused on the behaviour of large Wilson loops and the string tension. In parallel with the use of the reduced model, we also used the traditional extrapolation methodology to obtain the same quantity. Our comparison was presented in Ref. [32]- [33], and showed that the extrapolation of the finite N string tension agreed within the percent level with the quantity obtained from the reduced model. The result also matched reasonably well with other determinations of the large N string tension [34]- [35]. Obtention of the continuum string tension is an elaborate procedure involving noise reduction techniques on the raw data, extrapolation to large size Wilson loops and then to the continuum limit. All of these steps, as well as the traditional methodology, involve extrapolations. This is always a dangerous arena, although it could hardly be a coincidence that the results match so well. Our initial focus on continuum results is understandable, since after all the lattice is primarily a method to obtain properties of the continuum formulated Yang-Mills theory. However, the reduction idea should be also valid for the lattice model itself, even away from criticality. In this paper we want precisely to look into that. We will focus upon simple observables of the lattice theory. Some are robust, highly precise, and devoid of any technical issues which might shade the validity of the result. Thus, rather than testing reduction by checking the preservation of the symmetry needed to validate the original proof, we will produce a direct test by comparing the volume independence of the corresponding observables. Furthermore, this is a quantitative comparison which can put limits on the errors committed. Ultimately, these errors are crucial for a precision determination of observables.
What do we know about the finite N corrections to reduction? Perturbation theory gives some information in this respect. The propagators of the TEK model turn out to be identical to those of a lattice of size ( √ N ) 4 [15]. This is reasonable since we are essentially mimicking the space-time degrees of freedom by those of the group. This seems more efficient than the equivalent counting for the QEK model, which suggests an effective volume of (N 1/4 ) 4 . Furthermore, if we opt for partial reduction in the twisted case, the effective length of the box is given by L √ N . In addition to the effect on the propagators, the Feynman rules for the vertices also adopt a peculiar form. The structure constants turn into complex phases which depend on the "effective momentum" degrees of freedom. It is precisely this structure which suppresses the non-planar diagrams of the theory. The choice of flux enters the explicit form of the phase factors through a particular rational number in the exponent. Since all the factors cancel out for planar diagrams we were guided into thinking that the choice of fluxes was irrelevant. Now we know that this choice might be crucial in guaranteeing the preservation of centre symmetry. Recently [36][37][38][39] we examined the problem with more detail for the simpler case of the 2+1 dimensional gauge theory. Physical observables seem to depend smoothly on the rational factor. This suggests a stronger form of volume independence valid at finite N . According to the previous considerations, as long as we preserve the effective spatial size (LN in 2 dimensions) and the rational factor in the exponent of the vertices, we can trade spatial degrees of freedom by group degrees of freedom after a suitable tuning of the fluxes. It would be interesting to examine the situation for the 3+1 dimensional theory.
An important new ingredient appeared on the scene after the scientific community generated interest into quantum field theories in non-commutative space-times (for a review see Ref. [40]). It turns out that the afore-mentioned Feynman rules of reduced models are a discretized version of those appearing in those type of theories, where the rational number mentioned earlier is just proportional to the non-commutativity parameter. Indeed, the continuous version of the rules were actually anticipated by a generalisation of the reduced model obtained by one of us [41], and which captures the non-commutativity in its formulation. All these new connections suggest that the results obtained at finite N might have an interpretation and a usefulness in their own right.
In the previous paragraphs we have set up the situation that acts as background of the present paper. In the next section we will present all the methodological details leading to our results. This includes the description of the models and the parameters involved in them, the observables that we are going to study and the details about the simulations. The following section contains the analysis of our results for Wilson loops. The paper ends with a short concluding section.
2 Pure gauge theory on a finite lattice

The Models
As mentioned in the introduction we will focus upon pure gauge lattice theory at finite volume. The dynamical variables of our theory are the link variables U µ (n) which are elements of the SU(N) group. The index µ ranges over the 4 directions of space-time 0, 1, 2, 3, while n runs through all points of an L 4 hypercubic lattice L. The partition function of the model is given by where the link variables are integrated with the Haar measure of the SU(N) group, and S is action of the model. In this paper we will choose the simplest version of the action, proposed by Wilson, where U µν (n) is the unitary matrix corresponding to an elementary plaquette starting at point n and living in the µ − ν plane. Its expression in terms of the link variables is whereμ is the unit vector in the µ direction. The lattice coupling b is the lattice counterpart of the inverse of 't Hooft coupling λ = g 2 N . Thus, the large N limit has to be taken keeping b fixed. In the previous formulas, we assumed periodic boundary conditions for the variables: U µ (n + Lm) = U µ (n), for any integer vector m. Thus, the lattice is rather a toroidal one. Wilson action is invariant under gauge transformations of the link variables: where Ω(n) are arbitrary SU(N) matrices. All physical observables are required to be gauge invariant. This property allows to relax the periodicity condition to U µ (n + Lm) = Ω(n, m)U µ (n)Ω † (n +μ, m) (2.5) To make this formula well-defined, the same gauge transformation has to result from two different ways to reach the same replica point. Since matrices do not necessarily commute, this implies consistency conditions to be satisfied by the gauge transformation matrices Ω(n, m). In particular, one can consider n being the origin and m =μ +ν, and construct the gauge transformation in terms of the gauge transformations for m =μ and m =ν. The consistency condition reads: with z µν = exp{2πin µν /N } an element of the centre of the SU(N) group. Notice that if n µν = 0 for all planes, one can choose all Ω(n,μ) to be the identity and we recover periodic boundary conditions. On the other hand, if some n µν = 0 we get the so-called twisted boundary conditions introduced by 't Hooft [42]. The integer n µν = −n νµ can be interpreted as a discrete flux modulo N through the µ − ν plane. Now one can perform a change of variables to new link variables which are periodic at the expense of introducing a phase factor at certain plaquettes. The Wilson action now adopts the form S W (U µ ) = −bN n∈L µ =ν z µν (n)Tr(U µν (n)) (2.7) The z µν (n) factors can be moved around depending on the change of variables, but their product over each plane cannot be changed. One possibility is to make all factors equal to one except for a single plaquette sitting at the corner of each plane. It is quite clear that for large volumes this has little influence in the local dynamics, but at small volumes it has an important effect. An extreme situation occurs in reduced models. The lattice consists of a single point L = 1. One can drop the n dependence of links and plaquettes and one is left with a system involving only d SU(N) matrices U µ . If one takes periodic boundary conditions one has the EK model. If, on the contrary, one chooses twisted boundary conditions, one has the TEK model. The latter has more freedom because one can take different values for the fluxes n µν having different properties. For a class of values there exist configurations having zero-action. Thus, in these cases, there exist matrices U µ = Γ µ such that they satisfy These are the twist-eaters [43]- [44]. We will also impose the condition that the Γ µ matrices generate an irreducible algebra. The interested reader can consult the literature to see the conditions imposed on the values of n µν to achieve this goal (see for example Ref. [45]). Rather than working with the most general case we will concentrate on a rather simple situation which tries to preserve as much symmetry as possible among the different directions: the symmetric twist. The situation occurs whenever N is the square of an integer N =L 2 . Since, our goal is the large N limit, this restriction is not problematic. Furthermore, we will usually takeL to be a prime number to avoid Z(L) from having proper subgroups. In summary, for the symmetric twist case the twist factors for each plane are: whereL = √ N and k is an integer defined moduloL. Choosing k = 0 we recover the case of periodic boundary conditions.
For completeness let us specify the form of the Γ µ matrices for a symmetric twist with arbitrary k = 0. To do so, we first introduce twoL ×L matrices PL and QL whose non-zero matrix elements are given by They are SU(L) matrices satisfying The four N × N matrices Γ µ with N =L 2 , are then given by the direct product of PL and QL as follows with IL theL ×L unit matrix.

Observables
Now let us consider the main observables. The most natural gauge invariant observables are Wilson loops. A path on the lattice is a finite sequence of oriented links, such that the endpoint of one element in the sequence coincides with the origin of the next element. For each path we can construct a unitary matrix as a product of the associated link variables following, left to right, the order of the sequence. The reverse path is associated with the inverse matrix. A closed path is a path such that the endpoint of the last element in the sequence coincides with the origin of the first element. The trace of the corresponding unitary matrix is just the corresponding Wilson loop , which is gauge invariant. It is clear that the simplest non-trivial closed path is an elementary plaquette P µν (n). The associated unitary matrix U µν (n) ≡ U (P µν (n)) has the expression given earlier (Eq. 2.3). The corresponding expectation value E is the most precise lattice observable where the averaging is over positions, orientations and configurations. In pure gauge theories with Wilson action, the quantity is connected with the mean value of the action as follows where V is the lattice volume and d the space-time dimension.
For the case of twisted boundary conditions the change of variables into periodic link variables transforms the plaquette observable as follows: where the centre element z µν (n) is the same one appearing in the action. The above replacement in the formula for E gives the correct expression with twisted boundary conditions. A priori, E is a function of the different parameters entering the simulation E(b, N, L, n µν ). Given its local character this quantity is expected to have a well-defined infinite volume limit which is independent of the boundary conditions: E(b, N ). On general grounds one expects it to have also a well-defined large N limit E ∞ (b), with corrections that go as powers of 1/N 2 .
Apart from the plaquette, one can consider other Wilson loops. Here we will focus on those associated to a rectangular path of size R × T . The corresponding expectation value will be named W (R, T ). Once more, the quantity has to be modified in the standard variables for twisted boundary conditions. The modification is just to multiply the expression by the product of all z µν (n) factors for the plaquettes contained in the rectangle.
Just as for the plaquette, the expectation value will depend on b, N , L and n µν , but also depends on R and T . One expects a well-defined thermodynamic limit as L goes to ∞ at fixed R and T . For obvious reasons, the finite volume corrections should be larger as R and T get closer to L. These corrections should depend on the boundary conditions. An extreme situation occurs for the reduced model. In that case, the loop extends beyond the size of the box: Nevertheless, the statement of volume independence at large N implies that even loops that extend beyond the size of the box should behave as those obtained at infinite volume.

Methodology
The numerical results presented in this paper were obtained using a heat-bath Monte Carlo simulation. For the Wilson action case, each link variable U µ (n) has a distribution determined by the part of the action containing it, which has the form where V µ (n) is constructed in terms of the neighbouring links. Thus, at each local update one should generate a new link according to the previous distribution. However, the constraint that U µ (n) belongs to SU(N) complicates this task. As has been demonstrated in [46], a SU(N) link variable can be updated by successive multiplication of SU(2) submatrices. However, the standard Creutz' s heat bath-algorithm has a very low acceptance rate for SU(N) gauge theories with N larger than 8. To avoid this, we have used a heat-bath algorithm proposed by Fabricius and Haan [47] that significantly improves the acceptance rate for larger values of N. Apart from this, there are two technical points which significantly speed up the SU(2) heat bath manipulations.

Suppose
A is an SU(2) sub-matrix, and we want to change the link variable as follows U µ (n) ′ = AU µ (n). During the updating step, we should also store W = U µ (n)V µ (n), namely W ′ = AW , which only requires O(N ) arithmetic.
2. An SU(N) matrix has N (N − 1)/2 SU (2) sub-matrices (i,j). For odd(even) N, there are n 1 = N (N − 1) sets having n 2 = (N − 1)/2 (N/2) (i,j) pairs, which can be manipulated simultaneously. For example, for N=5, we have the following 5 sets of two pairs: If the computer supports vector or multi-core coding, we should make full use of these n 2 independent manipulations.
For the case of the reduced model (having L = 1), the previous methodology fails, because the distribution is not given by Eq. 2.18. The part of the action containing U µ (n) depends quadratically rather than linearly on it. However, one can return to a linear dependence by introducing a Gaussian random matrix, which after integrating over it reproduces the original form [47]. Thus, we can alternate the updating of the new random matrix with the previous heat-bath update.
Finally a link update is obtained by applying the SU(2) update for each of the N (N − 1)/2 SU(2) subgroups in SU(N). A sweep is obtained by applying an update for all the links on the lattice, followed by 5 over-relaxation updates.
To reduce autocorrelations we typically measured configurations separated by 100 sweeps. In any case, the final errors were estimated by a jack-knife method with much larger separation among groups. Another important point concerns thermalization. Typically we discarded the initial 10% of our configurations. We monitored several quantities to ensure that we observed no transient effects in the resulting data. Related to this point there is the choice of initial configuration. For the TEK model we always initiated thermalization from the cold configuration that minimises the action at weak coupling U µ = Γ µ . We observed no sign of phase transition behaviour except for low values of b and small N . Indeed, some of our results extend below the transition point between strong and weak coupling, which was estimated to be b = 0.3596 (2) in Ref. [48]. Results below this value of b lie then in a metastable phase (at infinite N ). Nonetheless, it seems that as N grows the probability of tunnelling decreases considerably and there are no signs of phase flips even down to b = 0.35. Whether the tests of volume independence are performed in a metastable region or not is of no particular concern, as was the case for the partial reduction studies of Ref. [29]. The same applies for the other phases presented in the same reference. In all our simulations we have monitored the behaviour of open loops which are not invariant under the Z( √ N ) 4 symmetry of the reduced model and found them to be compatible with zero within errors. As mentioned in the introduction, the choice of flux k is obviously crucial for this result. Our results, therefore, provide a test that the criteria, presented in Ref. [31] to avoid symmetry breaking, works well.
Finally, a comment about the total statistics and computer resources employed in our results. The number of different simulation parameters involved in our work are in the hundreds, and have been accumulated over several years. Typically, the number of configurations used for our large volume (L=16,32) periodic boundary conditions results is the range 300-600, and have been generated with INSAM clusters at Hiroshima University. For the TEK model we have used typically 6000 configurations at each b value at N=841, and 2000 configurations at N=1369. The most extensive and computer-wise demanding results were obtained with the Hitachi SR16000 supercomputers at KEK and YITP. Some small scale results were obtained with the clusters and servers at IFT.

General considerations about volume independence
Expectation values of the lattice observables depend on the conditions of the simulations, namely the values of b, L, N and the twist tensor n µν . Restricting ourselves to the symmetric twists we may write O(b, N, L, k) for a generic observable O. In the thermodynamic or infinite volume limit, local observables should have a well defined value: Notice that in that limit the dependence on the boundary conditions drops out. In general, we expect the infinite volume observable to have a well defined large N limit: The question now poses itself about the commutativity of the order of both limits. What is the result if we take large N limit first, and then the infinite volume limit? The concept of volume independence [10] gives a striking answer to this question. Taken at face value it predicts lim Even before it was formulated this hypothesis, as such, was known to be wrong. Its validity depends upon the value of b and the dimensionality of space-time. While presumably correct in two dimensions it is certainly not true in higher dimensions and sufficiently large values of b. Still the hypothesis is correct for b in the strong coupling region b < b c (N ) [10] [49]. The point brought into the discussion by Ref. [31] is about the dependence of the statement upon the value of k. No doubt that the choice of k has a crucial influence on the results for small L. The proposal of Ref. [15] of choosing k coprime with √ N avoids the problems found for k = 0 (periodic boundary conditions) at small enough couplings (large b).
A modified hypothesis for the periodic boundary condition (k = 0) case, known as partial reduction, was introduced by Narayanan and Neuberger [29]. In view of their results, they conjectured that volume independence remains true provided L > L c (b), with L c (b) a growing and asymptotically divergent function of b. This idea can be combined with the use of twisted boundary conditions. Not surprisingly, the effective value of L c (b), as obtained in finite N studies, seems to depend also on k [30].
However, numerical studies [17]- [18] disproved the initial hypothesis that it is enough to take k coprime with √ N =L to restore volume independence, as formulated in Eq. 3.3, at all values of the coupling. In view of this, a modified volume independence proposal was put forward in Ref. [31]: This implies that the flux has to be scaled asL is sent to infinity. According to the proposal, based on centre symmetry breaking arguments, the precise value of kL is not relevant provided kL/L andkL/L remain always larger than a certain threshold (∼ 0.1).
Herek is an integer such that kk = 1 modL. The restrictions of working with symmetric twists and values of N that are the square of an integer are presumably not indispensable. They emerge from the simplicity of maintaining as much isotropy as possible among the different directions of space-time.
The validity of any volume independence hypothesis has always been tested indirectly by verifying the Z 4 (N ) symmetry condition assumed in Eguchi and Kawai proof. Although this is valid road-map, nowadays it is possible to explore directly how a given particular observable behaves as a function of the different quantities involved. This has the advantage of informing us of the rate at which the asymptotic limits are attained, which is a crucial piece of information in rendering volume independence a useful tool in numerical studies. In what follows we will present the results of our analysis taking as observables the best measured lattice observables: the plaquette and other small Wilson loops. One of the good things of the reduction or volume-independence hypothesis is that it is valid for the lattice theory, not only the continuum limit. If reduction takes place in the scaling region, this will be inherited by the continuum limit.  E(b, N, L, k). The infinite volume coefficients up to order n = 3 have been computed by several authors [50]- [52]. The corresponding values in 4 dimensions arê

Behaviour at very weak coupling
The finite volume corrections for periodic boundary conditions started to be studied long time ago. The main difficulty is the presence of infinitely many gauge inequivalent configurations with zero-action: the torons [13]. If one expands around the trivial classical vacuum A µ = 0 (which dominates the path integral for d=4 and large N ), one must separate the fluctuations into those with zero momentum and the rest. The non-zero momentum contributions to the plaquette were studied in Ref. [53]- [54]. They adopt the form: The corresponding contribution of the zero-momentum degrees of freedom was computed in Ref. [55]. In 4 dimensions it becomes: It is clear that the leading order finite volume corrections do not cancel each other and do not go to zero as N goes to infinity: This kills the volume independence hypothesis for periodic boundary conditions at weak coupling, where the calculation is reliable. For twisted boundary conditions the calculation is very different since the action has isolated minimum action solutions. In this case the result is given by [13][55] Hence, there is no volume dependence at this order. This result, as well as that of periodic boundary conditions, can be easily derived since the expectation value of the plaquette can be obtained by differentiating the partition function with respect to b. To this order what matters is just the counting of gaussian and quartic fluctuations.
In the same references the volume dependence of rectangular R × T Wilson loops to leading order in perturbation theory was also studied. The infinite volume coefficients up to order λ 2 were computed in Ref. [51]. For loop sizes much smaller than the lattice size, the leading finite volume correction from non-zero momentum is while the zero momentum one is Again they do not cancel each other, leaving a finite volume correction at all N . On the contrary the result for twisted boundary conditions takes the form The first term is just the result of periodic boundary conditions at N = ∞ and a lattice volume of N 2 L 4 . Thus as N grows we approach the infinite volume result irrespective of the value of L. The second term goes like 1/N 2 in the large N limit, but plays an important role in reducing finite volume corrections. This comes because the leading L −4 corrections cancel each other between the first and second terms, so that the leading finite volume dependence in the sum goes like N −2 L −6 . Notice that, due to the same reason, adding the p = 0 contribution to the k = 0 coefficients of the right-hand side does not alter the result. In summary, large N volume independence holds at this order.
In Ref. [15] the present authors studied the structure of the Feynman rules for the TEK model (L = 1) and argued that volume independence should hold at all orders of perturbation theory. This result trivially extends to large N , k = 0 models for any L. The rules are essentially a lattice version of those of non-commutative field theory. The finite N corrections appear in two ways. On one side the propagators are just the lattice propagators in a box of size L √ N . This relates finite N effects with finite volume effects. However, in addition the finite N effect comes in as a non-planar diagram contribution, which is generically exponentially suppressed with N , but might induce also power-like suppressions. Furthermore, the coefficients depend on the choice of the flux k. To quantify finite N and finite volume corrections, including the effect of k, a numerical analysis of the O(λ 2 ) coefficients would be quite interesting. This study is currently under way [56].
In the absence of explicit perturbative calculations, we might use numerical simulations at small coupling in order to quantify the size of finite volume corrections at weak coupling (large values of b). In this work we used limited resources, by taking L = 1, 2, 4 and N = 49 for both periodic and twisted boundary conditions. A more precise study will be performed to test the higher order analytic calculations.
Hence, we will compare our numerical results with the predictions of perturbation theory at infinite volume. The results for the plaquette expectation value at b ≥ 2 and k = 0 are quite compatible with the perturbative formulas. For example, the perturbative prediction at infinite volume at b = 2 and N = 49 is 0.936152, while the measured value for tions. We already found that there are sizeable finite volume corrections of order λ = 1/b. However, even if we subtract out this contribution finite size effects are still considerably larger than for TEK. In Fig. 1 we show the case of L = 2 k = 0, and also the L = 1 k = 0, which has been divided by 16 to fit it into the same plot. It is also very interesting to look at what happens for larger loops. In contrast with the case of the plaquette, the result at lowest order has a non-vanishing finite size correction which, as we saw earlier, is of order 1/N 2 . The question is whether there is numerical evidence that higher orders violate this rule and produce large corrections. Again we have looked at the TEK model at N = 49 and various b ≥ 2 values. We measured the mean values of the Wilson loops and subtracted from it the perturbative formula for the infinite volume Wilson loop, up and including order λ 2 . Hence, the difference δW (R, T ) measures the leading corrections to volume independence. Asymptotically this must be dominated by the term linear in 1/b which can be computed analytically, and which for very large N and R, T ≪ √ N is given by Thus, the correction for R = T grows with the fourth power of (R/ √ N ). As expected √ N acts as the effective size of the box. Thus, in Fig.2 we plot our result for δW (R, R)/R 4 at N = 49 with twisted boundary conditions and k = 3. The straight lines are the leading order perturbative contribution. Since, R = T = 3, 4 are not much smaller thanL = 7 the slopes for the various R change by a factor of 2. The actual data lie always below the corresponding straight line, implying that higher powers of 1/b produce a dependence of the opposite sign. Indeed, one can easily describe the data points by simply adding a 1/b 2 correction with a coefficient determined by a fit. The corresponding line is also plotted as the blue curve for R = T = 4.  Fig. 3a as a function of 1/N 2 . The dependence on N is clearly visible, being much larger than the individual errors, which are too small to be seen in the plot. A good fit is obtained by a second degree polynomial in 1/N 2 . The χ 2 per degree of freedom is 0.248, for 9 points and 3 parameters. The coefficient of the linear term is 0.960(3). The constant coefficient in the fit provides the extrapolated value of the plaquette at N = ∞, and its value is E ∞ (0.36) = 0.558012 (12). To test the L dependence we repeated the analysis in an 8 4 and 4 4 periodic box. The results of the 8 4 box are numerically close to those of 16 4 , but the difference is much larger than the errors and clearly seen in the plot. For N = 8 the difference is 3.41(70) 10 −4 and for N = 16 it is 2.42(50) 10 −4 . A similar three parameter fit describes the data very well and gives an extrapolation of 0.558114(67). Indeed, even fixing the extrapolated value to be equal to the L = 16 one, we get a fit with a chi square per degree of freedom which is smaller than 1. Thus, our results are consistent with the basic claim that the volume dependence drops to zero in the large N limit.
In contrast we also show the results of the L = 4 periodic lattice. For this smaller volume we have explored values of N = 8, 10,12,14,16,25,49,81. A good polynomial fit can be obtained which undoubtedly extrapolates to a very different value in the large N limit. Thus, volume independence fails in this case. This is precisely what we expected on the basis of the results of Narayanan and Neuberger [29]. At b = 0.36 volume independence should only work for L ≥ L c (0.36) ∼ 8. The important role of boundary conditions, as claimed in Refs. [14]- [15], is certainly supported by our data. The black point in Fig. 3a is our result with the TEK model, namely a one point box (L = 1) with symmetric twisted boundary conditions, flux k = 11 and group rank N = 1369. From a single simulation we get E ∞ (0.36) = 0.558019 (11), which is in perfect agreement with the extrapolated result. Obviously, it is to be expected that the same applies if one uses partially reduced situations with L > 1, provided one uses appropriate boundary conditions. Indeed, we studied the L = 2 and L = 4 cases with symmetric twisted boundary conditions and N = 9, 16, 25, 36, 49, 81 (with k values of 1, 1, 2, 1, 2, 2 respectively). Here the N dependence of the result is quite obvious, but follows the same curve as for the L = 16 periodic boundary conditions one. We take this curve to give the N dependence at infinite volume E (0.36, N ). The leading correction, as mentioned earlier, is a 1/N 2 term with a coefficient close to one. We use this fact to display the data in a different fashion, which allows to get a better grasp of the precision involved. Plotting E − 1/N 2 the errors of the L = 16, 8 data are clearly seen in the plot (Fig. 3b). Notice the small scale on the y axis, showing the high precision of the data. The curves shown are the afore-mentioned fits to the periodic boundary condition data in the range N ∈ [8,16]. Once again we see that the curve goes through the TEK point at This shows that L = 4 is close enough to infinite volume when twisted boundary conditions are used even for N as small as 9. Our L = 2 data with tbc do not disrupt the picture. Indeed, a simultaneous fit to the L = 2, 4 twisted and L = 16 periodic box has chi square per degree of freedom smaller than 1.
The previous data provide the most precise test of volume independence to date. Errors in the plaquette are of order 10 −5 , and to that level the Twisted Eguchi-Kawai model matches the extrapolated result. At this level of precision the 1/N 2 correction is still sizeable up to N ∼ 300. We emphasise that this agreement is putting to test the validity of reduction in the non-perturbative regime. The value of b = 0.36 lies in the region in which scaling is best observed and corresponds to a lattice spacing of 0.1 in Λ MS units. Thus, it is hard to argue against the fact that nonperturbative phenomena contribute corrections to the plaquette higher than 10 −5 .
In any case, as argued in the introduction, our goal goes beyond the yes/no answer to volume reduction. We want to investigate the size and nature of the errors affecting physical quantities at large but finite N . For the case of the plaquette we have seen that there is a leading 1/N 2 correction which is pretty much volume independent provided one uses sufficient large sizes and/or appropriate boundary conditions. For example, at L = 16 and N = 16 the plaquette is E(0.36, 16) = 0.561700 (6), which is still relatively far from the infinite N value of 0.558. Using twisted boundary conditions (k = 1) and L = 4 the result is E(0.36, 16) = 0.561712 (20). For L = 2 one has E(0.36, 16) = 0.561748(84). Thus the three results are perfectly consistent with each other. Unfortunately the L = 1 twisted result for this N cannot be obtained due to the vicinity to the strong coupling phase. The plaquette shows frequent jumps to a low value around 0.4. As mentioned in the previous section these flips are suppressed at larger values of N .
We have seen that the N -dependence of the plaquette expectation value for large lattice sizes matches with the result of small L = 2, 4 sizes and symmetric twists. In principle, there is no known reason why the 1/N 2 corrections should not depend on L and the value of k. For the case of the TEK model (L = 1) it is important to know how big these corrections are, in order to optimise the values of N used to reproduce infinite volume results within a given precision. To investigate this matter, we studied the TEK model for N = 81, 121, 289, 529, 841, 1369 corresponding toL = 9, 11, 17, 23, 29, 37. The values of the flux adopted were k = 2, 3, 5, 7, 9, 11 respectively. These choices are dictated by our criteria to avoid symmetry breaking. The resulting values of the plaquette are collected in Fig. 4. We plot E(0.36, N ) − 0.96/N 2 according to the value of the leading correction obtained from the L = 16 periodic boundary condition data. This correction only affects sizably the results of N = 81 and 121. What we see is that all the results are consistent with each other within errors. That shows the approximate universality of the 1/N 2 correction. On the negative side of the x-axis at arbitrary locations we placed the prediction for the infinite volume, infinite N , plaquette E ∞ (0.36) obtained by extrapolation of the L = 16 and L = 8 periodic boundary conditions data. The results are consistent, and give a best estimate E ∞ (0.36) = 0.558002 (5). The different values and errors are also displayed in Table 1. We should emphasise that the choice of k is important. If we approach those values for which symmetry breaking is observed, significant departures are observed. For example, taking N = 81 and k = 1, the plaquette becomes 0.5554 (1). Although this is less than 1% away from the result at k = 2, it is many standard deviations away given the precision of our data.

Results for other small loops and for other values of b
Having seen that reduction works at the level of the tiny statistical errors for the plaquette, we should also explore other quantities. In principle, extended quantities should be more sensitive to finite volume corrections than the plaquette. Hence, we explored also the results for square loops of linear size R (The R = 1 case is just the plaquette). The relative errors grow considerably with R. Typically, errors remain more or less constant, but the value of the observable drops considerably; at R = 4 by a factor of a hundred. Since, we want to test the observables with the minimum amount of manipulation we restrict our analysis to R ≤ 4. The methodology is similar to the one for the plaquette. The L = 16 results are fitted to a quadratic polynomial in 1/N 2 allowing extrapolation to large N. The coefficient of the 1/N 2 is 1.04(1), 0.44(2) and 0.14(1) for R = 2, 3, 4 respectively. For the larger loops the quartic coefficient is compatible with zero. The large N extrapolated values are collected in Table 1 together with the plaquette values. The fits are good and match with the result of the TEK model. As an example, the situation for the largest loop W (4, 4) is shown in Fig. 5. Notice that, as a consequence of the small expectation value of the loop, the 1/N 2 correction gives a quite significant contribution in relative terms (of order 10% at N = 16). The results for L = 8 periodic boundary conditions give compatible extrapolations for all loops except for R = 4 as seen clearly from the figure.
Hence, the evidence that at b = 0.36 reduction is at work with errors of order 10 −5 extends to Wilson loops of size up to 4 × 4. A less extensive study has been carried at other values of b. For example, at b = 0.37 and L = 16 we also measured the same observables in the range N ∈ [8,16]. Again the data is beautifully described by a second degree polynomial in 1/N 2 . The linear coefficient is now 0.789 (5). This gives an estimate of E ∞ (0.37) = 0.578978 (17), which compares quite well with the direct measurements of the TEK model at N = 1369 and N = 841 which were E ∞ (0.37) = 0.578961 (6) and E ∞ (0.37) = 0.578954 (7) respectively. The situation is depicted in Fig. 6 in which we To the level of the errors it is clear that the 1/N 2 correction is rather large. We also analysed the characteristic 1/N 2 dependence for the TEK model (L = 1). Obviously for that one has to simulate small values of N , since for N > 300 the correction is smaller than the statistical errors. Furthermore, it is to be expected that the coefficient depends on the flux k. Thus, we simulated the model for various (k, N =L 2 ) combinations. In all cases we computed the coefficient C 1 given by The results are plotted in Fig. 7 as a function ofk/L (we recall thatkk = 1 modL). Within the allowed region (> 0.1) the coefficient lies between −2 and 2, implying that there are no unexpectedly large 1/N 2 corrections in any case. Furthermore, all values havingk L > 0.25 are consistent with the infinite volume 1/N 2 coefficient marked as a horizontal line in the plot. For bigger loops the 1/N 2 is rather a finite volume correction, which should grow as the loop size approaches the effective lattice sizeL. Nonetheless, the equivalent C 1 coefficient remains within reasonable limits. For example, for 4 × 4 loops the coefficient remains in the band [−6, 10] fork/L > 0.1.
The most important qualitative difference of the b = 0.37 data with respect to the b = 0.36 one is that the L = 8 periodic boundary condition data does not seem to extrapolate to the same value. The result is 0.57935(7) which is not terribly far away, but inconsistent within errors. Indeed, this was to be expected on the basis of the Narayanan and Neuberger results, since at b = 0.37 the L = 8 size lies outside the unbroken symmetry region.
The previous procedure was repeated for various other values of b. Two estimates of the large N infinite volume Wilson loop expectation values can be given. One is the large N extrapolation of the results on an 16 4 periodic box. The other comes from simulating the TEK model at various large values of N . The TEK results are consistent among themselves and compatible with the ones coming from large N extrapolation. Furthermore, they tend to have smaller errors. Table 2 summarises the results.

Global Fit to small loop expectation values
In the previous subsections we have analysed the validity of the reduction idea by comparing the results of infinite volume gauge theories with those obtained for the TEK model. Our first tests were obtained in a region dominated by perturbation theory, and then we moved to the opposite edge of the weak coupling region. It is tempting to use all the results to  If we analyse all plaquette expectation values obtained for the TEK model, it is remarkable that all our results for b ≤ 0.385 fall neatly into a straight line in the variable 1/b I , where b I = b * E ∞ (b). The improved coupling b I was introduced by Parisi [57] and used extensively in Ref. [34]. The function 0.89661(32) − 0.006803(7)/b I approaches the numerical values at the level of 0.0001. Nonetheless, the fit cannot be considered good, since the errors in the plaquette values are one order of magnitude smaller.
This led us to look for another parameterization, which could not only fit the measured values but also match with the perturbative expression up to three loops studied earlier. It was not hard to find a solution as a ratio of two polynomials of third degree in 1/b I . The constant term of both polynomials is 1. Of the remaining six parameters, three are fixed by the perturbative formula. Fitting the remaining three parameters to the TEK values we get a fit with chi-square per degree of freedom equal to 0.22. Given the nice properties of the formula and the high precision of the data, this seems quite remarkable.
We tried to extend the previous formula to the determination of E(b, N ). As mentioned earlier, in addition to the L = 16 data used for the extrapolations presented in previous subsections, we do have a good amount of data coming from our determination of the large N string tension [32]- [33]. Altogether, we have 88 points for L = 16 and L = 32 with periodic boundary conditions, various values of b and N ranging from 3 to 16. A good fit is obtained to all of them by replacing the 3 parameters of the previous parameterization by a second degree polynomial in 1/N 2 . The constant term is fixed to the parameters of the large N fit. Altogether, there are 6 new parameters, although some are clearly redundant and one can obtain a good fit having a chi square per degree of freedom of ∼ 0.7 with only 4 parameters. Furthermore, 1/N 4 terms are only necessary if we include the results at N = 3 and N = 4.
Summarising, we may writeĒ The denominator has exactly the same form with different coefficients a ′ nm which are completely determined by the a nm and the coefficients of the perturbative expansion of the plaquette to order 1/b 3 (we leave the details to the reader). Notice that Eq. 3.18 gives E(b, N ) in an implicit way, since b I = bĒ(b, N ). Nevertheless, the formula allows one to solve analytically forĒ(b, N ) as a function on b with a formula that involves radicals.
In Fig. 8 we show the difference of our plaquette measurements and the function Eq The conclusion of this result is that the whole data set is perfectly consistent. All values can be fitted with a simple function whose large N behaviour has been determined entirely in terms of the TEK measured values. The finite N values are dominated by the 1/N 2 correction which describes the data from a relatively small value of N on.
Furthermore, having this formula allows one to estimate the plaquette expectation value at intermediate values. For example, we could compare our formula with the measurements of other authors at large volumes and different couplings. We always found agreement to the level of ±0.00002. Nonetheless, one must be cautious about using our formula for extrapolation outside of the fitting range. In principle, for large b this is less worrisome, since the function incorporates the perturbative behaviour up to order 1/b 3 . To give a bird's eye view of our formula and the data we show the plaquette value E as a function of b for several values of N in Fig. 9. To the scale of the figure our best fit curves go right through the centres. For comparison we added also the curve reconstructed from the 34 perturbative coefficients obtained recently for the SU(3) case [58]. At this scale this curve also fits the points and you only start to see visual deviations of both curves for low b values.
The same procedure can be extended to square loops up to 4 × 4. There are a bunch of small differences in our procedure and also in the quality of the results. We considered more reasonable to fit the logarithm of the Wilson loops. This time, however, we lack the O(1/b 3 ) perturbative coefficients for all N , and we used only the results of Ref. [51]. In order to deal with this lack of information without enlarging the number of parameters we -25 - used a (3,2) Padé approximant with coefficients that are polynomials in 1/N 2 . The lowest order 3 parameters are determined in terms of the TEK data only and give good fits. The higher order terms in 1/N 2 are determined from the SU(N) data with poorer fits having chi-squares of 2 times the number of degrees of freedom (points-parameters). Given the small errors there are still tiny differences between the measured values and the fits. In Fig. 10 we give as an example the bird's eye view of the 4 × 4 Wilson loop, which at this scale shows no visible deviation. As explained earlier, the interest of the fit is that we have been able to make use of all the data obtained during the last few years by our group both for the TEK reduced model as for ordinary SU(N) lattice field theory at various values of N . The main conclusion is that everything fits into a picture in which the TEK results match nicely as the large N extrapolation of ordinary gauge theory at infinite volume for the wide range of couplings explored.

Conclusions
In this paper we have presented the results of an extensive analysis of the expectation values of small Wilson loops for SU(N) Yang-Mills theory on the lattice with Wilson action. Our analysis allows to test the dependence of these observables on the rank of the group N , the size of the box L and the boundary conditions. Our results provide direct evidence that,  under certain conditions, the large N results are independent of the lattice size L. For this to happen the choice of boundary conditions is crucial. Twisted boundary conditions are most effective in reducing the volume dependence at large N , as advocated in Ref. [14]- [15]. Our data are consistent with the claim of Ref. [31], that even the one point model (TEK) captures the infinite volume large N results if the large N limit is taken with a suitable choice of the fluxes. Hence, our conclusion is that an appropriate version of volume independence holds. Although this is not a mathematical proof, the observables studied can be determined with a high accuracy. Hence, the validity of the statement has been tested down to the level of the statistical errors of our data, which are of order 10 −5 .
We have analysed the region of large values of b = 1/λ where perturbation theory is a good approximation, but most of our results were obtained in the interval b ∈ [0.35, 0.385] typically used to extract continuum limit results. In any case, there is no evidence for any pathological behaviour at any value of b. Indeed, we were able to find a parameterization of the infinite volume plaquette and small square loop expectation values which encodes the known perturbative behaviour at weak coupling and fits all our measured values. Furthermore, the large N limit of this parameterization matches nicely with the results of the TEK reduced model. A large amount of data obtained over the last few years has gone into this analysis. As stressed in the introduction the goal was not only to check the volume independence hypothesis but also to estimate its corrections. At infinite volume the Wilson loops that we studied approach the large N limit with a leading correction that goes as 1/N 2 with a coefficient (depending on the value of b) of order 1. The correction is not small at the scale of the errors (∼ 10 −5 ) and only becomes comparable for values of N ∼ 200 − 300. Thus, performing simulations at various N and extrapolating can hardly be avoided. Furthermore, for the plaquette, the typical 1/N 2 corrections of the reduced model have similar size to those of SU(N) lattice gauge theory. Hence, a single simulation of the reduced model gets the correct result with orders of magnitude less degrees of freedom. In all our previous statements care has to be taken on the choice made for the flux integer k. However, our results do not show problematic dependencies provided one stays within the safe region proposed in Ref. [31].
For larger loops one must take into account the connection between the √ N and an effective lattice size which follows from perturbation theory. Hence, deviations are expected if the loop sizes become close to √ N . Here, we have shown that provided one remains away from this situation, things look pretty much like for the plaquette.
In this work we intended to verify the commutativity of the large N and large volume limits for typical and precise lattice quantities. Our test does not rely on centre symmetry and Eguchi-Kawai proof. Things are nevertheless consistent since our order parameters for the symmetry are compatible with zero in all our simulations. If reduction holds for the lattice model, this property should be inherited by the corresponding continuum limits. Indeed, this track was followed earlier [33] when we attacked the determination of the continuum string tension at large N . Although our results were strikingly consistent with volume independence, they employed sophisticated methods of noise reduction for large loops as well as extrapolation to the continuum limit. This paid a price from the point of view of simplicity, as well as considerably enlarging the errors. Hence, the present work can be seen as a necessary complement.