Large $N$ scaling and factorization in SU($N$) Yang-Mills gauge theory

The large $N$ limit of SU($N$) gauge theories is well understood in perturbation theory. Also non-perturbative lattice studies have yielded important positive evidence that 't Hooft's predictions are valid. We go far beyond the statistical and systematic precision of previous studies by making use of the Yang-Mills gradient flow and detailed Monte Carlo simulations of SU($N$) pure gauge theories in 4 dimensions. With results for $N=3,4,5,6,8$ we study the limit and the approach to it. We pay particular attention to observables which test the expected factorization in the large $N$ limit. The investigations are carried out both in the continuum limit and at finite lattice spacing. Large $N$ scaling is verified non-perturbatively and with high precision; in particular, factorization is confirmed. For quantities which only probe distances below the typical confinement length scale, the coefficients of the $1/N$ expansion are of $\mathrm{O}(1)$, but we found that large (smoothed) Wilson loops have rather large $\mathrm{O}(1/N^2)$ corrections. The exact size of such corrections does, of course, also depend on what is kept fixed when the limit is taken.


Introduction
An interesting approach to study quantum chromodynamics (QCD) is to consider the order of the gauge group N as a free parameter. As shown by 't Hooft [1], by taking the large N limit of the perturbative weak coupling expansion, the theory simplifies in many ways, and in fact one can treat theories at finite N as corrections in the "small" parameter 1/N . Moreover, the large N expansion predicts that quark loop effects are suppressed by a power of 1/N , so that the weak coupling expansion of large N QCD is dominated by planar diagrams with purely gluonic internal loops. All of these rather remarkable properties of large N QCD, make it an interesting theory to study not only from the theoretical perspective, but also from a practical point of view, as results for real world QCD could be obtained by considering corrections to the N = ∞ theory which are parametrized by powers of 1/N . Although this 1/N scaling is obtained perturbatively, lattice computations provide evidence that it also holds at the non-perturbative level, both in D = 4 spacetime dimensions [2][3][4][5][6][7][8] and in D = 3 [9][10][11][12][13]. The evidence is usually based on complicated observables, where typically one needs to project onto ground states by large Euclidean times. It is then difficult to obtain high precision at various N in order to verify 't Hooft scaling with good confidence. Let us stress the fact that the validity of the 1/N scaling, beyond the weak coupling expansion, is not a trivial statement. Hence, it is desirable to test it by means of lattice simulations and with statistically and systematically very precise observables.
Perturbatively, if one carries on with the 't Hooft 1/N topological expansion, another simplification arises, which has to do with the property of factorization where the O i are local gauge invariant or Wilson loop operators, and the leading correction scales as 1/N 2 in the pure gauge theory, which we focus on for the rest of this work. Eq. (1.1) has several consequences, as it tells us that in the large N limit, the dominant part of a correlator is the disconnected one. In particular, when O 1 = O 2 , this means that fluctuations are suppressed; and as discussed in Ref. [14], this fact can be put in analogy with the classical limit of a quantum theory, where 1/N plays the role of . Related to this is also the concept of the "master field", i.e., the idea that the path integral is dominated by a single gauge configuration (or rather a gauge orbit) [15,16]. Although these ideas triggered hope to find the solution of large N QCD, such an analytical solution is still lacking today. The situation in the Yang-Mills theories is similar in this respect to two-dimensional SU(N )×SU(N ) spin models [17], while for O(N ) models and CP(N ) models the large N limit is solvable and one can therefore really carry out the expansion [18][19][20]. One more aspect where Eq. (1.1) plays a crucial role has to do with the idea of volume independence, which starting from the work of the authors in Ref. [21], has been used in the lattice formulation to study the large N limit of the Yang-Mills theory by performing simulations in small spacetime volumes [22], and even in single site lattices, provided a clever choice of boundary conditions [23][24][25] is made. Clearly, the possibility to compute observables on a single site lattice makes simulations of SU(N ) Yang-Mills theory at large N more accessible, as there is a significant compensation of the extra cost for increasing N by the much smaller number of lattice sites.
The above indicates that factorization is not only relevant in the theoretical context, but also on the practical level, as it is a requirement for the single site lattice simulations to be valid. To be more precise, the equivalence between the single site and the infinite volume theory is argued for on the basis of the Makenko-Migdal loop equations [26] on the lattice. As originally shown in Ref. [21], the loop equations in both theories are equivalent, provided that the product of the expectation value of the Wilson loops factorize as stated in Eq. (1.1). Additionally, we would like to point out that important physics is contained in the corrections to factorization. The most obvious one is that glueball masses are obtained from the connected correlation functions of Wilson loops.
The previous discussion motivates the search for a non-perturbative proof, beyond the realm of weak coupling perturbation theory. Several authors have investigated factorization beyond perturbation theory [27][28][29][30], but no verification has been carried through using the non-perturbative framework provided by numerical simulations of lattice gauge theories. We here fill this gap. In addition, the observables that we consider make it possible to study the large N scaling up to very high precision, and hence address the important issue of the size of the corrections to N = ∞.
This paper is organized as follows, in Sec. 2 we present the observables that are used both to check the large N scaling, as well as factorization. In Sec. 3 we discuss different ways of defining the large N limit and in particular the two choices we made for our investigation. In Sec. 4 we describe the ensembles and lattice parameters used for the simulations and in Sec. 5 we present our results, both at finite lattice spacing, and in the continuum limit. We finish with a short summary of the results.

Observables
The basic observables we consider are the Yang-Mills action density E(t) at positive flow time [31] (defined below) and rectangular Wilson loop operators where C is a closed rectangular path in space-time, and P denotes the path ordering operator. The normalization factor 1/N is included in the definition of W C in order to have a finite large N limit -already at tree-level. Wilson loops have singularities which have to be removed before the continuum limit can be taken. In particular, for our square Wilson loops, one must remove not just the "perimeter" divergences but also "corner" divergences [32][33][34]. One way to proceed is to consider Creutz ratios [35], which however, for loops of large size in lattice units, suffer from small signal to noise ratios. As this would compromise our desire for a precision test, we work instead with smooth Wilson loops. The smoothing is provided by the Yang-Mills gradient flow [36,37]. It evolves the gauge fields A µ (x) according to the flow equation where the dimension two parameter t is known as the flow time. The loops at positive flow time are then simply Choosing 8t to be of a typical QCD size, say of the order of the square of the inverse string tension, they benefit from small statistical errors even for large loops [38]. In particular, their variance remains finite in the continuum limit. That property is a particular manifestation of the most important feature of observables which are built from the smoothed gauge fields B µ (t, x): they are renormalized operators at positive flow time t [31,39]. In other words, there is no renormalization scheme or scale dependence beyond t and the continuum limit is unambiguous and well defined. Even the action density is finite. It will be one of our observables.

The gradient flow coupling at large N
The gradient flow can also be used to define a renormalized coupling [37]. Using the perturbative expansion of the Yang-Mills energy density at positive flow time E(t) , one has that 3 γ E + 52 9 − 3 ln 3 is N independent. With this definition, we can then define a scale by setting the renormalized couplingλ GF to a given value. A convenient choice for SU(3) is the reference scale t 0 [37], which corresponds to a value of the coupling such that t 2 E(t) | t=t 0 = 0.3. This particular choice can be generalized to SU(N ) if the right hand side is modified so that it has the correct scaling with N .
Clearly, we also want the definition to remain what it is for N = 3. Thus, Eq. (2.5), suggests to define t 0 implicitly by the equation [2] t 2 E(t) | t=t 0 = 0.1125 for all N .

Smooth Wilson loops
The favourable properties of smooth Wilson loops have already been exploited in the literature, as for example to estimate the string tension at small values of t in Refs. [6,39], or to study the large N phase transition in the eigenvalue spectrum of the Wilson loop matrices [40]. For our purpose, the limit of small t is not required, as the smooth loops are used to test factorization and the large N limit for well defined renormalized observables, regardless of their relation to the operators at t = 0.
In the end, we study the large N limit of square Wilson loops, i.e. for loops where the path C in Eq. (2.3) is given by a square of size R × R. In order to take the large N limit, the loops are matched at different N relating their size to the scale t 0 introduced in the previous section. More precisely, the large N and continuum limits are taken for loops of size R c = √ 8ct 0 (see Figure 1), where the smoothing parameter t = ct 0 , and c is a constant parameter.
To be more precise, let us denote a square loop with one of its corners at the spacetime point (x 0 , x) and extending only in space as W (t, x 0 , x, R). Its expectation value is independent of the position x due to translation invariance and only depends on the parameter c. In our notation we separate time and space-coordinates, as we will later use lattices with different boundary conditions in time (open b.c.) and space (periodic b.c.). While the independence on x is exact, x 0 has to be sufficiently far away from the boundary for Eq. (2.7) to hold. Similarly, we define which corresponds to the expectation value of the product of a Wilson loop with itself.

Observables to test factorization
In order to investigate the property of factorization from Eq. (1.1), we define several observables based on the Yang-Mills action density and the smooth Wilson loops at positive flow time. They are constructed such that factorization implies that they vanish as N → ∞. First we consider the simplest case of the observable G W defined in terms of the smooth Wilson loops as Then, we consider observables built from the space integral of the smooth Wilson loops and the Yang-Mills action density. We define 1 with the factor 1/t 3/2 0 rendering H O (c) dimensionless, and where O is either a smooth Wilson loop, or the Yang-Mills action density. Notice that H is a type of susceptibility, as we are integrating over the contributions from the correlation function of O at different distances. The integration does not extend over x 0 due to our choice of boundary conditions and x 0 is again supposed to be far away from the time-boundaries. In comparison to the simple observable G W , this probes longer distances, but introduces also more noise and affects the statistical errors in the measurements. Nonetheless, as will be shown in Sec. 5, the statistical precision that can be achieved for H O remains good. In particular, we will consider H E (c), defined by inserting We remind the reader of our choice (2.13)

Finite volume
For a numerical test, we need to choose a finite volume. We chose our parameters such that L/ √ 8t 0 ≈ 3.3. Table 1 shows the actual values used in our simulations. Since L is thus approximately constant, it is omitted as an argument of the observables. We note that the large N limit and factorization can be tested in infinite or in finite volume. To be on the safe side, we chose the latter, even though we are not far from the infinite volume limit for most observables.

Defining the approach to the large N limit
The complete definition of a quantum field theory involves a regularization (here Wilson's lattice theory) as well as a non-trivial renormalisation before the regulator can be removed. Although this is usually not discussed, quantitative statements about the approach to the large N limit, such as the ones we are seeking here, do depend on the renormalisation scheme if the renormalisation scheme defines which quantity is held fixed as we take N → ∞.
While the O(1/N 2 ) corrections depend on these details, the true limit is expected to be unique in the following sense. It is independent of the scheme, as long as the 't Hooft couplingλ s (µ) = Nḡ 2 s (µ) in any scheme is kept fixed as one takes the limit. This statement becomes most transparent when we replace couplings by the associated Λ-parameters, Now any renormalization group invariant quantity O of mass dimension n, has a large N limit and Examples for n = 1 are glueball masses and t 0 defined above is a RGI scale with n = −2. When the observable O depends on external momenta or coordinates, they have to be fixed in units of Λ in a specified scheme, e.g. Λ MS , when taking N → ∞.
Due to the existence of the limit eq. (3.2), we may also scale distances with respect to any one particular reference scale (choice of O). In our numerical work we have chosen t 0 , eq. (2.6), because of its high precision.
The preceding discussion is about the continuum theory. It thus saliently assumes that first we take the continuum limit at finite N and then we perform N → ∞. However, we may also proceed in the opposite order: first take the large N limit at fixed lattice spacing and then send the lattice spacing to zero. 2 Let us briefly discuss that this order of limits is indeed the same as above; the limits are interchangeable.

Large N limit at fixed lattice spacing
The existence of the large N limit at fixed finite lattice spacing is expected due to the following consideration. We start from the Lambda-parameter, Λ lat in the lattice minimal subtraction scheme, which satisfies eq. (3.1) with µ = 1/a and λ lat (µ) = λ 0 in terms of the lattice spacing, a, and the bare coupling, λ 0 = N g 2 0 . In fact, having a specific scheme, the lat-scheme, we can give the more detailed formula, . Eq. (3.6) shows that the large N limit can be taken at fixed bare coupling which is equivalent to fixed lattice spacing a. Apart from O(1/N 2 ) terms in c 1 and higher order terms, fixed lattice spacing is the same as fixed Λ lat and therefore also fixed Λ s in other schemes. See also an early discussion of Λ MS /Λ lat including its N dependence [42]. In general, taking the large N limit at fixed lattice spacing has to be followed by the a → 0 limit at N = ∞. However, when we investigate factorization, the second step is not expected to be necessary. This is because the perturbative proof of factorization holds in the lattice regularization [43] at finite a. If factorization holds non-perturbatively we thus also expect eq. (2.13) at any fixed a. In any case, verifying eq. (2.13) at arbitrary finite lattice spacing implies that it holds in the continuum limit.
Note also that even the large N limit of divergent quantities, such as Wilson loops at t = 0, is expected to exist. A high precision numerical test has recently been performed [8].

Lattice details
In this section we give the details of our lattice simulations. We simulate the pure gauge theory with N = 3, 4, 5, 6, 8 at several lattice spacings. The lattice action is the Wilson gauge action and we use open boundary conditions in the time direction [44]. The simulations are performed using a combination of heatbath and overrelaxation local updates using the Cabibbo-Marinari strategy [45] to refresh the SU(N ) matrices. The ratio of overrelaxation to heatbath updates is fixed to L/(2a).
For convenience, we present the values of the lattice spacing, as well as lattice sizes in physical units by assigning a value to t 0 such that √ t 0 = 0.166 fm. This choice is motivated by the result in SU (3) for √ 8t 0 /r 0 = 0.941 (7) [46] and the value of the reference scale r 0 ≈ 0.5 fm [47]. Notice that this choice is somewhat arbitrary, as apart from the missing quark loops, for N = 3 the theory cannot be directly identified with Nature.
The parameters of the simulations are displayed in Table 1. The configurations used for the measurements are a subset of those reported in Ref. [2] for all ensembles except for those at N = 3, 8, and for the finest lattice spacings in the case of N = 4, 5.
As announced above, all the lattices considered in Table 1 are of approximately the same spatial size L ≈ 1.55 fm. In addition, we have used two additional ensembles with L ≈ 2.35 fm at the coarsest lattice spacing (a ≈ 0.1 fm) for N = 4, 5 in order to check for effects due to small variations in the volume. Notice that for the ensembles which have been reported in Ref. [2], we have a very large number of measurements for the Yang-Mills action density.
The flow equations are integrated using a third order Runge-Kutta integrator [37] and the observables are measured at intervals ∆t of t of ∆t/a 2 ≈ 2−3×10 −2 . Afterwards, they are interpolated using a second order polynomial in order to obtain their values at arbitrary t. The action density is defined exactly as in [37], using the clover discretization and it is measured from t = 0 up to t ≈ 1.2 t 0 . The loops, W (c), are measured only in the vicinity of t = ct 0 , with c = 1/2, 1, 9/4, and then interpolated to the exact value of t. For the loops, one has to do an additional interpolation to R c , and since their statistical precision is very high, one has to be careful with small potential systematic effects. The details of this interpolation were already presented in Ref. [49].
We end this section with the precise definition of the observables introduced in Eqs. (2.7)-(2.10) on the lattice. First, for the Wilson loops we use translation invariance in the form   (N ) we give the inverse lattice coupling β = 2N 2 /λ 0 , the dimensions of the lattice, the approximate lattice spacing using √ t 0 = 0.166 fm followed by the number N W meas of measurements used for the computation of the smooth Wilson loops, and N E meas for the action density, eq. (2.4). In the second to last column we present the values of t 0 /a 2 : * taken from Ref. [48] and * * taken from Ref. [2].
where · corresponds to the estimator of the true expectation value computed on the lattice.
In order to compute H W , we define so that and we proceed in a similar way to define H E after replacing W (ct 0 , x 0 , x, R c ) by t 2 E(ct 0 , x 0 , x)| t=ct 0 . The parameter d is introduced to deal with the systematic effects from the open boundary conditions. It is chosen in a similar way as described in Ref. [50], so that the effects coming from the boundaries are negligible with respect to the statistical error in the bulk.

Large N scaling
In order to test and provide a precise verification of 1/N 2 scaling, we analyse our results for W (c) and for the gradient flow couplingλ GF . Let us first discuss our results for the latter.

The gradient flow coupling
In Figure 2 we plotλ GF as a function of t for several gauge groups and different lattice spacings. Within the scale of the plot, the results are hard to distinguish for all gauge groups, which already shows the small size of the N dependent corrections. While at t = t 0 indepencence of N is enforced by Eq. 2.5 or equivalentlȳ the different N curves remain remarkably close when t is a factor of 5 away from t 0 . At a closer look, corrections to N = ∞ are present and the data agrees very well with a polynomial in 1/N 2 as expected. We verified this by interpolating the data to several values of t in a regular interval from t = 0.1 t 0 to t = 1.1 t 0 , and then taking the large N limit once at a fixed lattice spacing and once in the continuum. As can be observed in Figure 2, cut-off effects are large at small t. At t = 0.1 t 0 the relative difference between the results at the finest lattice spacing (a ≈ 0.05 fm) and at the coarsest (a ≈ 0.1 fm) one, is around 20%; while the errors in the measurements themselves is at the per-mill level. The situation is better at larger values of t, so let us first focus on values of t/t 0 ≥ 0.3, where the relative size of cut-off effects is reduced tenfold, when compared to the case at t/t 0 = 0.1. In Figure 3 we show a plot of the continuum extrapolation ofλ GF at t/t 0 = 0.8 and the large N extrapolation both at finite lattice spacing and in the continuum. In order to be able to use the dataset at N = 8, in addition to the continuum limit extrapolations, we consider a 2 /t 0 = 0.2091, the value on ensemble A(8) 2 . We then interpolated the results for all the other gauge groups to that lattice resolution.
On the left panel of Figure 3 we show the continuum limit extrapolations. The strategy chosen for the extrapolation is the following: all continuum extrapolations are performed by linear fits in a 2 /t 0 to those data which satisfy a 2 /t 0 ≤ 1/4 (default fit). Such a restriction has been well motivated in Ref. [51] for N = 3 and we find smaller discretization effects for larger N . As an estimate of the systematic uncertainty associated with this choice, we perform a second fit linear in a 2 with a data point at larger a 2 ; if the latter fit does not have a good χ 2 we add an a 4 /t 2 term to the fit-function (control fit). If necessary, the error of the default fit is enlarged until it covers the full 1-σ band of the control fit. The uncertainties of the continuum limit points are usually dominated by the systematics which arises from different fits and which is not necessarily independent for different N .  Figure 2:λ GF as a function of t for several values of N and a (see Table 1). In the lower plot, we present a closer look at the small t region.
All values of χ 2 /dof are excellent except for SU(4), where we obtain a value of 2.2 and 2.7 for the linear and quadratic fits respectively. After performing the fits in a 2 /t 0 , on the right panel of Figure 3 we plot the large N extrapolations both in the continuum and at finite lattice spacing. As discussed, N = 8 is available only at finite lattice spacing, where in addition, errors are much smaller due to the fact that we performed an interpolation instead of the continuum limit extrapolation. The large N extrapolation uses the form As seen in Figure 3, the fit to the function Y is excellent, with a χ 2 /dof = 1.02 at finite lattice spacing; for the continuum points we do not consider χ 2 since the errors are strongly correlated due to the dominating systematic uncertainty of the continuum extrapolations. Notice also that the results suggest that cut-off effects decrease with increasing N .
As an example of results at smaller t, we show our analysis at t/t 0 = 0.4 in Figure 4. In this case, the magnitude of the cut-off effects is larger, but the same analysis as before can be carried out.
As mentioned earlier, dealing withλ GF at values of t/t 0 ≤ 0.3 presents a bigger challenge, so one cannot reach the same level of accuracy as the results presented in this section. However, we have proceeded to do a similar analysis for such small values of t, including corrections of higher order in a 2 /t 0 . Details are found in Appendix B.
From the above analysis, we find that the large N dependence ofλ GF is in excellent agreement with the 1/N 2 scaling predicted by the 't Hooft perturbative expansion. Moreover, defining we can determine the "distance" between SU(3) and SU(∞). In the continuum, at t/t 0 = 0.8 and t/t 0 = 0.4 we find η(1/9) = 1.1% and 2.8% respectively. Note also that the large N limit is taken at fixed t 0 and therefore η ≡ 0 at t = t 0 by definition.
To account for this effect, we also fit , and define δ in a similar way to η. The results, together with those obtained for η are displayed in Table 2. Let us remark that the individual errors in our measurements are below the per-mill level, so we can confidently quantify these percent level deviations between SU(3) and SU(∞).
The magnitude of the 1/N 2 corrections can be read off from the coefficients a 1 and a 2 collected in Table 2, together with those of the smooth Wilson loops which we discuss next.

Smooth Wilson loops
We have determined the smooth Wilson loops at three different values of c, i.e. c = 1/2, 1, 9/4. As in the case ofλ GF , we are interested in the large N scaling at finite lattice spacing and in the continuum. The strategy for the continuum limit fits is the same as forλ GF . The fits for the loops at different c are qualitatively similar, so in Figure 5 we show the results at c = 1 only.     Once again, to quantify the magnitude of the finite N corrections, we collect in Table 2 the values of a 1 and a 2 from the fit to Y . We observe that the relative magnitude of them grow at larger values of c (or t equivalently). Similarly, the deviation between SU(3) and SU(∞) also grows up to a value of η(1/9) = 0.1 when c = 9/4. In all cases we find an excellent fit to Y (the values of χ 2 /dof are reported in Table 2).

Factorization
In order to verify the property of factorization from Eq. (1.1), we take the large N limit of the observables defined in Sec. 2.3. The large N limits are taken in a similar way as described earlier, but we modify the parametrization of the large N fitting function for convenience, so that For the continuum limit extrapolations we use the same strategy as for W and forλ GF , and in all cases, the data can be fitted very well with a linear or quadratic polynomial in a 2 /t 0 . We also check for effects caused by variations of L/ √ 8t in all observables. As discussed in Appendix A, we find that H W at c = 1/2 and c = 1, are potentially affected by large effects. We tried to include them as a systematic error on the measurements, but this yields errors which are too large to be of interest as a test of factorization. Hence, we present only results for H W at c = 9/4 . Let us first discuss our results for H E (1). On the right panel of Figure 6 we show the large N fits both in the continuum and at a finite lattice spacing. The fits are excellent, which provides yet again confirmation of the scaling in powers of 1/N 2 . It is worth mentioning that at finite lattice spacing, where results are very precise, we find that a quadratic fit in 1/N 2 , excluding the SU(3) point, extrapolates to  Table 3. At finite lattice spacing we include also the parameters from the fit excluding SU(3). Concerning the N → ∞ limit itself, the extrapolated value is within two standard deviations from zero in the worst case. Notice that at finite lattice spacing, the errors in the extrapolation are two orders of magnitude smaller than the value of H E (1) at N = 3. To further validate factorization, an additional fit is performed for which b 0 = 0 is fixed, and only b 1 and b 2 are fitted to the data. This enforces factorization, so the value of χ 2 /dof from the fit can be used to asses the validity of the assumption (L * , C * in Table 3). To summarize, for H E (1) we find excellent agreement with factorization in the continuum, and a deviation compatible with two standard deviations in the worst case at finite lattice spacing, still statistically consistent with factorization.
We now turn to the smooth Wilson loops. In Figure 7 we display the results of the continuum and large N fits for G W (1). The parameters of the extrapolations at the three values of c are displayed in Table 3. Also in this case, we find that SU (3) can be used as a validation point, and if it is excluded from the fit, it agrees with the extrapolating function within two standard deviations at c = 1, and within one standard deviation at the remaining values of c.
For the fits with factorization enforced by fixing b 0 = 0, the values of χ 2 /dof are also excellent. These values, together with those of b 0 reported in Table 3, give us confidence on the validity of factorization. Notice that the errors at large N are at least one order of magnitude smaller than the value at SU(3) itself. Concerning the finite N corrections, comparing the loops at different values of c, we observe that those at large c are characterized by large coefficients in front of the 1/N 2 and 1/N 4 correction terms.
Yet another interesting question is whether loops at fixed t but different R have obs.    (3), and C: in the continuum. Additionally, we fit the data to a function with b 0 = 0, so that factorization is imposed at finite lattice spacing L* and in the continuum C*. In this case, the value of χ 2 /dof validates this hypothesis. different finite N corrections. We explore this issue at the end of the section. Let us first look at H W in Figure 8. The N → ∞ limits are less than two standard deviations away from zero. Inspecting the continuum limit fits, we observe that had we taken the three-point extrapolation for N = 6 as our central result, the central value would have been close to the upper end of the error bar in Figure 8 and the 1/N 2 extrapolation in full agreement with factorization. In other words, one should not look too much at the central value but at the full range of the error, as always.

Loop size dependence
Finally, let us explore how the finite N corrections to factorization change when the size of the loop is increased at a fixed value of the smoothing parameter t. For a given value of t, we consider square loops of size R(ξ) = ξ √ 8t. Given the finite size    of the lattices, we use the loops measured at c = 1/2 (t = t 0 /2), so that we can consider larger values of ξ. Thus, at a fixed value of t = t 0 /2, we defineŴ andĜ W in a similar way as W and G W , but in this case, as a function of ξ instead of c. In addition to the already presented results at ξ = 1, we also measuredŴ andĜ W at ξ = 1.25, 1.5, 1.75 and 2. The coefficients obtained for the large N fits at finite lattice spacing are displayed in Table 4 as a function of ξ. We observe that in the case of the loops themselves, the coefficients of the 1/N 2 expansion do not change significantly with ξ, while those ofĜ W grow rapidly for larger loops. In fact, they grow exponentially fast as shown in Figure 9. At finite N , larger loops are much further away from N → ∞ than smaller loops.

Conclusions
We have taken the large N limit of a few observable of SU(N ) pure gauge theories numerically defining all dimensionfull quantities in units of t 0 . This means that we held t 0 , or equivalently the couplingλ GF , eq. (2.5), at a low energy fixed in defining the approach to the limit. As explained in Sect. 3, the precise magnitude of 1/N 2 corrections do depend on this choice. For each quantity, the continuum limit was taken before the large N limit, but we have also investigated large N scaling at finite lattice spacing, defined by a 2 /t 0 =constant.
In both cases we find that finite N observables are very well and very precisely described by a leading order term and corrections ∼ 1/N 2 and ∼ 1/N 4 . We recall for example Figure 4 where the excellent precision, in particular at finite a, is visible. In the same way, factorization has been confirmed very precisely. Of course, a numerical computation cannot substitute a mathematical proof, but our results make it very implausible that anything goes wrong with the large N limit in general, or factorization in particular.
However, the magnitude of corrections to the large N limit is more complex. We found a strong dependence on the physical size of the observables. For example, we considered R × R Wilson loops smoothed with a smoothing radius of size again √ 8t = R. Table 2 shows the deviation, η, of SU(3) from SU(∞) of these smooth loops to increase from 3% at a loop-size of r = 0.2 fm to 10% at R = 1 fm.
When we increase the loop size R at fixed smoothing radius √ 8t = 0.23 fm from R = 0.23 fm to R = 0.5 fm, the corrections η(1/9) ≈ a 1 /9 (with a 1 from Table  4 or Figure 9) grow from 4% to more than 30%. The growth with R of the finite N corrections to factorization is even more dramatic as seen on the right panel in Figure 9. These large corrections may also contribute to the fact that one has to go to very large N to approach the large N limit in the 1-point model [8,52]. Of course, the dominating effect is expected to be that the color degrees of freedom provide the effective size of the system in that model.
In summary, large N scaling is confirmed with high precision, but corrections to large distance observables can be substantial. One thus has to be careful when deriving quantitative information from large N considerations in gauge theories. are within the statistical uncertainty. To showcase this, in Figure 10 we display a plot of H W (c) for c = 1/2 and c = 9/4. Clearly, at the smaller c, volume effects are much larger. Both at c = 1/2 and c = 1 we find that including the volume effects as a systematic correction is difficult with just two points in L. Attempting it in a conservative manner produces errors which are too large to check for the large N scaling at a similar precision as for the rest of observables, so we include only H W (9/4) in the analysis. We note that we do not have an explanation why small c appears to be more difficult than a large one.
Bλ GF at small t As stressed in Sec. 5.1.1, the continuum extrapolations become more difficult at smaller values of t. Due to the large cut-off effects, we find that a linear extrapolation in a 2 /t 0 does not parametrize the data adequately, even using only the finest lattice points. For that reason, we include in the fits the O(a 4 /t 2 0 ) and the O(a 6 /t 3 0 ) corrections. The fit strategy is similar to the one used for the rest observables, except that higher degree polynomials are used. Briefly, the central point is obtained by performing a quadratic fit in a 2 /t 0 using those data for which a 2 /t 0 ≤ 1/4. Then, a second fit is performed, either using a quadratic function, or a cubic one if the value of χ 2 /dof of the first one is too large. Finally, the error is chosen so that it covers the full 1-σ band of both fits. We show our results at t/t 0 = 0.2 in Figure 11. As expected, the errors in our continuum extrapolations are larger than those obtained at larger values of t, but within the errors, the large N extrapolation is perfectly consistent with a polynomial in 1/N 2 . At finite lattice spacing, the large N extrapolation is cleaner, and once again it shows and excellent agreement with the 't Hooft expansion.