Numerical stochastic perturbation theory applied to the twisted Eguchi-Kawai model

We present the results of an exploratory study of the numerical stochastic perturbation theory (NSPT) applied to the four dimensional twisted Eguchi-Kawai (TEK) model. We employ a Kramers type algorithm based on the Generalized Hybrid Molecular Dynamics (GHMD) algorithm. We have computed the perturbative expansion of square Wilson loops up to $O(g^8)$. The results of the first two coefficients (up to $O(g^4)$) have a high precision and match well with the exact values. The next two coefficients can be determined and even extrapolated to large $N$, where they should coincide with the corresponding coefficients for ordinary Yang-Mills theory on an infinite lattice. Our analysis shows the behaviour of the probability distribution for each coefficient tending to Gaussian for larger $N$. The results allow us to establish the requirements to extend this analysis to much higher order.


Introduction
Gauge theories in the limit of infinite number of colours (N −→ ∞) are very interesting theoretically [1]. They are simpler than their finite N counterparts, but share most of the fascinating properties of the latter. Nonetheless, their understanding remains a challenge. The relevance of this goal is also given by the fact that they seem a point of contact with other approaches [2]. A good deal of their difficulty lies in the non-perturbative character of most of its properties. The standard first principles approach to this kind of problems is the lattice formulation of quantum field theory [3]. However, in contrast with what happens in perturbation theory, the large N limit involves an extrapolation and seems harder than the finite N study. Only recently this line of approach has led to trustworthy computations (for a review see [4]), pioneered by the works of Teper and collaborators [5][6][7].
Fortunately, there is a very specific simplification which emerges when studying large N gauge theories on the lattice. This is the so-called Eguchi-Kawai (EK) reduction [8]. According to this result finite volume corrections are subleading in the large N limit. Hence, there is the possibility that the dynamics of the large N gauge theories is captured by a matrix model. There are several proposals that have been put forward to transform this possibility into a reality [9][10][11][12][13]. Here, we will be concerned with one of the early proposals called Twisted Eguchi-Kawai model (TEK for short) [14,15], and introduced by two of the present authors. Indeed, it was used shortly after EK proposal to compute the string tension at large N [16,17]. Lately this has led to a calculation of this quantity with at least compatible precision to other methods [18]. Several precise tests of the reduction mechanism have been obtained recently for other quantities [19], providing a strong verification of the validity and usefulness of this approach.
The previous paragraph justifies our interest in the TEK model. The finite N corrections are different for the matrix model and for the ordinary gauge theory, and their size and nature are very important from the practical point of view, in order to make this approach competitive computationally. Curiously, these corrections have a theoretically interesting interpretation in terms of the so-called non-commutative gauge theories [20]. Their Lagrangian and Feynman rules first appeared in the literature [21] when looking for a continuous generalization of the TEK model, and before they emerged from the mathematical construction of non-commutative geometries [22]. A more direct connection appears as a result of Morita duality on the non-commutative torus. This shows that ordinary gauge theories with twisted boundary conditions (TBC) a la 't Hooft [23] are particular cases of non-commutative field theories. The TEK model is nothing but the volume reduced version of a gauge theory with twisted boundary conditions on the lattice. TBC are characterized by a collection of integer-valued fluxes, and their appropriate choice has been found to have crucial impact upon the size of the corrections and the absence of phase transitions.
The last ingredient entering in this work is perturbation theory. Although our major interest when using the lattice approach was in studying the non-perturbative aspects of the theory, there is also interest in understanding the theory from the perturbative side. At the least it offers us a method to analytically determine certain observables and estimating the size and nature of its N dependence. In this spirit a recent perturbative calculation of Wilson loops in lattice gauge theories with twisted boundary conditions (including the TEK model) was addressed [24].
In addition, a new set of ideas has focused in explaining the usefulness of perturbation theory in understanding also non-perturbative contributions. This takes its most extreme form in the concept of resurgence (see for example [25,26] and references therein). At the least, as was already known, the large order behaviour of the perturbative coefficients is associated with non-perturbative aspects such as the action of other saddles. Very interesting results have been obtained in this spirit for large N matrix theories [26,27]. For example, it has been identified that the large N phase transition of the Gross-Witten-Wadia unitary matrix model [28,29] is governed by non-trivial saddles and the action can be reconstructed via the resurgence on the asymptotic expansion about the vacuum [30].
To close the circle, recently a new method has arisen that allows the numerical computation of high order coefficients in the perturbative expansion. It goes under the name of Numerical Stochastic perturbation theory (NSPT) and was pioneered by the works of Di Renzo and collaborators [31][32][33][34][35][36] (for earlier developments on stochastic quantization [37] and stochastic perturbation theory, see also [38][39][40]). Recently, this methodology has led to remarkable results about high order coefficients in SU(3) gauge theories [41][42][43][44][45][46][47][48][49]. As the memory size for fields expanded perturbatively is proportional to the highest perturbation order involved in the computation in NSPT, the authors of [45,47,48] have used twisted boundary conditions aware that they lead to reduced finite volume dependence.
The Langevin equation has been employed in the stochastic quantization, and NSPT is based on the perturbative expansion of this Langevin equation. For lattice theories various Monte Carlo algorithms, such as Kramers, hybrid molecular dynamics (HMD), and hybrid Monte Carlo (HMC) algorithms, have been developed and used. Recently extensive studies have been carried out for the NSPT algorithms [50,51]. Especially the HMD based algorithm for NSPT is easy to implement by perturbatively expanding the HMD/HMC codes used for non-perturbative simulations and one could benefit from various numerical integration schemes for the molecular dynamics part to reduce the finite step size error of the integration. The NSPT algorithm based on the HMD/HMC algorithm could open the way to efficiently estimating very high order coefficients to detect non-perturbative effects from perturbation theory numerically.
Our work is a first attempt to apply this methodology to matrix models. For the time being, our work is mostly exploratory but a necessary step before any attempt of a higher order and larger N study. Nonetheless, apart from showing that the method works with an incredibly high precision for the low-order coefficients, we also present results which extend to higher order the previous analytical results [24]. In addition, our work provides interesting information about the probability distribution of the perturbative coefficients. In particular, we have studied the dependence of the cumulants of these distributions with respect to the different parameters: the value of N and of the fluxes, the size of the loops and the order of perturbation theory. These results allow us to determine the necessary computational requirements for any further extension of these studies.
The layout of the paper is the following. In section 2, after introducing the TEK model, we explain the application of the NSPT algorithm based on the HMD algorithm to the TEK model. In Section 3 we present the numerical results for the perturbative coefficients of the Wilson loops. We employ (N, k) = (16, 1), (49,2), (121, 3) for SU(N ) and flux parameter k of the TEK model. The coefficients are computed up to four-loop (O(g 8 )) and the first two-loop coefficients are compared to the analytic values. The probability distribution of these estimates is investigated and used to explore the feasibility for extending our results to higher orders. We estimate the numerical computational cost of NSPT for the TEK model at large values of N and large perturbative orders in section 4. Possible improvements on the algorithm are also discussed. In the last section we summarize the main results of the paper. Two technical points are included in appendices.

TEK model and NSPT
In this section we start by briefly introducing the TEK model together with the gauge fixing method. Then we recall the Hybrid Molecular Dynamics (HMD) algorithm for nonperturbative simulations of the TEK model. This algorithm is then perturbatively expanded to derive the equation of motion for the NSPT algorithm.
The computational cost is estimated in terms of the highest order N trunc in the perturbative expansion involved in the algorithm and the matrix size N of SU(N ).

TEK model
The partition function of the TEK model in four-dimensions is defined by where U µ are SU(N ) matrices. We choose the symmetric twist characterized by N =L 2 and z µν given by where k is an integer coprime withL (= √ N ). The bare coupling constant g is defined through β = 1/g 2 .
The Wilson loop operator for an R × T rectangle in the µ-ν plane is defined by In this paper we restrict ourselves to square loops R = T . The analytic perturbative expansion of the observables proceeds by expanding the link matrices around the classical vacuum as where the matrices Γ µ define the classical vacuum of eq. (2.2) and have the following property A particular solution (up to multiplication by a phase) is singled out by an appropriately chosen gauge fixing condition, accompanied by the corresponding ghost term. NSPT also needs the gauge fixing via the stochastic gauge fixing method to stabilize the runaway trajectory of the field variables [52]. Here we define the gauge fixing functional for the Landau gauge condition from which we can extract the stochastic gauge fixing contribution in NSPT [31,32]. The gauge fixing functional F [G] is defined by where G ∈ SU(N ) is the gauge transformation matrix. The Landau gauge condition is achieved by maximizing the functional, by application of several iterations of the form: where α is a parameter. When this iteration process is transformed into a continuous evolution equation with a fictitious time t, we obtain the following differential equation: where w(t) is an auxiliary Hermitian-traceless matrix tracing the steepest descent trajectory. This evolution equation can be combined with the NSPT process to implement the stochastic gauge fixing.

NSPT for the TEK model
NSPT is based on stochastic quantization [37]. This amounts to writing down a stochastic differential equation of the Langevin type for the field variables of a target system, which relaxes to the equilibrium probability density at large times [38][39][40]. The NSPT method is obtained by expanding these field variables in powers of the coupling constant and casting the original Langevin equation into a tower of equations, one for each power. Observables that are analytic functions of the field variables acquire a corresponding perturbative expansion. Each term in the expansion becomes an stochastic variable whose mean value at large times gives us the corresponding coefficient that we are after.
In this paper we employ the HMD based NSPT, which has been introduced in refs. [50,51]. This is easy to implement through modifications of existing codes of nonperturbative HMD or HMC algorithms, and is preferable to systematically improve the molecular dynamics (MD) part, reducing the error arising from the finite MD step size.
To explain the method in our case, we will begin by reviewing the generalized HMD (GHMD) algorithm for nonperturbative simulations. Then we will apply the perturbative expansion to the algorithm, to derive the HMD based NSPT algorithm.
Step 1: Set initial momentum P 0,µ from the Gaussian distribution as where the matrices T a are SU(N ) generators in the fundamental representation normalized as Tr(T a T b ) = 1 2 δ ab .
Step 2: Solve the MD equation defined bẏ from the initial state (P 0,µ , U 0,µ ) for a fixed interval in fictitious time. The dot˙represents the left-wise derivative with respect to the fictitious time. These evolution equations preserve the energy value E = H[P, U ].
Step 3: Store U µ as the configuration and set U 0,µ = U µ , then return to Step 1.
The HMD partition function is introduced by adding the canonical momentum variable P µ conjugate to U µ to the original partition function eq. (2.1). The HMD algorithm proceeds as described in Algorithm 1, where the equations in Step 2 are discretized in time. This discretization violates the exact energy conservation, which distorts the distribution shape. In order to ensure that the discretized Markov chain leads to the correct probability distribution the Monte Carlo process must satisfy time reversibility and area preservation in the MD evolution, for which the leapfrog type scheme is normally used. For nonperturbative simulations, the Metropolis test can be inserted between Step 1 and 2 to compensate the violation of energy conservation, yielding the HMC algorithm.
The full momentum refreshment in Step 1 will cause random walking in phase space so that the autocorrelation times becomes longer for the ensemble. To relax random walking behaviour, the generalized HMC algorithm (GHMC) has been introduced in ref. [53]. For the GHMC algorithm the momentum is partially refreshed in Step 1 as instead of eq. (2.17), where c 1 is a mixing parameter. P µ in the right hand-side is the solution of Step 2. A momentum reflection step is added after the rejection of the Metropolis test in the GHMC algorithm, while this step is absent in the GHMD algorithm. We can further include the gauge fixing process eq. (2.13) into the MD evolution equation.
With the leapfrog algorithm and the partial momentum refreshment, the updating algorithm for one step having a small ∆t interval, is given by where (P 0,µ , U 0,µ , w 0 ) is the initial state and (P 1,µ , U 1,µ , w 1 ) is the final state. Eq. (2.25) corresponds to the leapfrog evolution for ∆t, for which various higher-order schemes are available. For simplicity we will explain the 2nd order leapfrog scheme only. The gauge fixing evolution (2.12)-(2.14) is interleaved into the MD evolution [54,55] and the transformation in eq. (2.26) does not affect the gauge invariant observables, however, we introduce this to explain the stochastic gauge fixing for NSPT [31,32,52]. Taking the ∆t → 0 limit with c 1 = e −γ∆t , as shown in ref. [51], this evolution reduces tȯ where D µ w is defined as and ζ µ is a random noise satisfying The terms with w act as the gauge damping force. Eq. (2.27) corresponds to the Kramers equation and it with γ → +∞ (c 1 = 0) to the Langevin equation with gauge damping force.
In the GHMD algorithm Step 2 is obtained by repeating eqs. (2.25) and (2.26) until the total evolution time becomes a fixed value, and then applying eq. (2.24) to implement Step 1. When the evolution time in Step 2 reduces to ∆t, the method becomes Langevin (c 1 = 0) or Kramers (c 1 = 0) algorithm.
Having explained the GHMD algorithm, we now describe the corresponding NSPT algorithm. This follows by expanding (P µ , U µ , w) and the MD equation (2.22) or eqs. (2.24)-(2.26) in a power series in g. We will now describe the perturbative expansion for eqs. (2.24)-(2.26), because it is sufficient to write the simulation program. In order to simplify notation, we first define the -product as the convolution product of two perturbative series, given by where A, B, C are matrices and A (k) , B (k) , C (k) are the coefficient matrices.
The perturbative expansion for (P µ , U µ , w) is defined by where U (0) µ = Γ µ is kept fixed in NSPT as it is the perturbative vacuum. Rescaling the fictitious time and the gauge fixing parameter as t = t/g and α = g 2 α, and substituting eqs. (2.33)-(2.35) into eqs. (2.24)-(2.26), we obtain for each perturbative order k = 1, 2, . . . . The details of the perturbative expansion of the matrix exponential exp µ and Θ (k) are extracted from eqs. (2.21)-(2.22), and (2.10)-(2.11), respectively: (2.43) In the limit ∆t → 0, the equation of motion should conserve the energy eq. (2.16) order by order in perturbation theory. This energy conservation can be monitored during the simulation. The Wilson loop operator eq. (2.4) is expanded similarly: where we used a shorthand notation for matrix power with the -product such as The expectation values at large times of our W (k) µν (R, T ) yield the coefficients of the perturbative expansion of Wilson loops. For the latter we take the notation given in Algorithm 2 NSPT algorithm based on GHMD. The MD time step size ∆t = 1/N MD , momentum mixing parameter γ, and gauge fixing parameter α are given.
0 ) out to SU(N ) group and algebra.
ref. [24] where the first two coefficients have been computed analytically. If we consider for example an R × T Wilson loop in the µ-ν plane, the relation is as followŝ Notice that, since the perturbative expansion is in powers of 't Hooft coupling λ = g 2 N , the expectation value of W (k) µν for odd values of k has to vanish. A technical point is the necessity to correct for deviations from the unitarity constraint induced by the numerical round-off error from the finite precision of computer arithmetic. While imposing the traceless-Hermitian character of (P (k) µ , w (k) ) is easily done, the conditions on {U (k) µ } following from the unitarity of the link matrices can be imposed by the matrix logarithm scheme for the reunitarization [36]. The details of the group projection is described in Appendix B.
Let us conclude by summarizing our NSPT algorithm. As explained earlier, the finite time step induces a distortion in the probability distribution, which in NSPT cannot be corrected by a Metropolis test. Hence, it is preferable to employ a higher order integration scheme in the MD evolution. Although, for simplicity, we used the simple (2nd order) leapfrog scheme to explain the HMD based NSPT algorithm, in practice we employed the 4th order leapfrog scheme for eq. (2.37), while the evolution of w is not changed. The properties of 2nd and 4th order Omelyan-Mryglod-Folk schemes [56,57], and 2nd order Runge-Kutta scheme in NSPT have been investigated in [50,51]. We summarize the NSPT algorithm used in this paper in algorithm 2. This corresponds the Kramers type NSPT algorithm (called KSPT in [50]). The specific parameters are fixed to γ = 0.5 for the momentum refreshing parameter c 1 = e −γ∆t and α = 2 for the gauge damping parameter. The trajectories are performed over a fixed time t = 1, with a perturbative reunitarization step at the end. The Wilson loop coefficients are computed at the end of every trajectory.

Computational estimates
The maximum power of g that is studied N trunc is limited by the total computer memory available. For the SU(N ) TEK model, the total memory requirement is of O(N trunc N 2 ). Most of the computational time is spent in matrix multiplication, whose cost is of O(N 3 ). The cost of the convolutional -product is of O(N 2 trunc ) with a naive implementation. The evaluation of the perturbative matrix function requires one more factor of N trunc , yielding the cost of O(N 3 trunc ) (see Appendix A). The cost of the reunitarization we implemented is of O(N 4 trunc ) (see Appendix B). The energy difference in a trajectory ∆H is non-zero for a finite MD step size ∆t. The relation between the energy conservation violation ∆H and ∆t depends on the MD integration scheme. We assume that algorithm 2 with 4th order leapfrog scheme yields ∆H (k) ∼ ∆t 4 at each perturbation order [50], where ∆H (k) is the energy conservation violation for the perturbative coefficient of the total energy in NSPT. Since the energy is proportional to N 2 , to achieve a constant energy conservation violation, the number of steps N MD , therefore, should scale as N MD = √ N . The total computational cost of the NSPT algorithm truncated at N trunc for the SU(N ) for the reunitarization part. Our estimates of total CPU time do not take into account autocorrelation times. For that we need experimental studies on the statistical properties of the Monte Carlo data, which depends on models and algorithms. As far as the parameters we investigated, no sizable autocorrelation and parameter dependence are observed and the autocorrelation length is of O(1) in units of the trajectory length t = 1.

Numerical results
In this section we show the results for the perturbative coefficients of the Wilson loops with NSPT, and investigate the statistical properties of the distribution. We have implemented the NSPT algorithm as explained in the previous section. We accumulate the statistics for the perturbative coefficients of Wilson loops up to order g 8 (N trunc = 8). The mean value gives us the corresponding perturbative quantity, but the variance and higher cumulants allow us to estimate the required statistics to achieve a given error in these coefficients. Fortunately, we have analytic results for the one and two-loop coefficients to test our results at these orders, but we can extend these results two more orders in powers of λ. Furthermore, we monitor the coefficients at odd-orders O(g 2 +1 ) which should be zero for Wilson loops. This provides a test of the cancellation of the non-loop effect in NSPT. Indeed, our results are consistent with zero within the two standard deviation. Therefore, we will concentrate in giving the results of even-order coefficients.
The parameters and statistics are shown in table 1. N MD is the number of MD steps for unit trajectory t = 1. Statistical errors are estimated with the jackknife method after binning in 1000 trajectory samples. We discard the first ≥ 1000 trajectories to account for thermalization. Since we employ the 4th order leapfrog scheme for the MD integrator, a ∆t 4 = 1/N 4 MD dependence is expected in observables [50]. We also study the limit of vanishing step size for some representative cases. Given the pilot nature of our study we have concentrated in studying cases which have been analyzed in detail in the analytic calculations of ref. [24]. Thus we concentrated in the three values of N = 16, 49, 121. The lowest values are not so interesting from the point of view of approximating large N Yang-Mills theory at infinite volume. Corrections are expected to be large. However, only for a value as low as N = 16 can one see clearly two effects of the twisted boundary conditions: the breakdown of CP and of cubic-rotational invariance. The first phenomenon manifests itself in non-vanishing imaginary parts of Wilson loops. We included this case to test this phenomenon in our NSPT determinations. On the other hand having N = 49 and 121 allows us to test the dependence on N of both the physical and computational parameters. We also tested several values of the flux parameter k. For sufficiently large values of k/L the dependence is quite small. For large values of N and small values of k the dependence can be quite strong even in the perturbative calculation. Non-perturbatively this restriction is even more important to preserve a remnant of center symmetry [58] that validates reduction in the limit N → ∞.

Comparison with analytic calculations
As mentioned earlier the presence of the twist breaks part of the symmetries of the standard lattice symmetries with periodic boundary conditions. In particular, part of the cubic group is broken. This translates into a dependence of the coefficients of the     spectively. The mean values are plotted as a function of 1/N 4 MD and extrapolated linearly to vanishing step size. The linear dependence in this variable is the expectation for our 4th order leapfrog scheme for the MD integrator, and our data are consistent with this expectation. The data for the S 1 and S 2 planes are plotted and extrapolated separately (purple dashed : S 1 , green dotted : S 2 ). The analytic results, depicted by horizontal lines in the same plots, only predict differences among planes forŴ 2 . At all orders, the cubic symmetry breaking in the real parts is found to be small and of the order of the errors of our calculation. As mentioned earlier the imaginary parts beyond the leading order are non-zero and different for the two families of planes. Our results reproduce both features and match nicely with the analytic results for the first two coefficients.
All our results for (N, k) = (16, 1) and (49, 2), extrapolated linearly in 1/N 4 MD to zero, are collected in Tables 2 and 3. We also tabulate the analytic values from ref. [24] and χ 2 /DoF from the extrapolating fit. The violation of CP and cubic invariance due to the twist is well seen for the N = 16 case at the two loop level and the differences between S 1 S1 S2  and S 2 averages are consistent with the analytic values. For the (N, k) = (49, 2) case, the results of our NSPT analysis are also consistent with analytic values. However, the violation of CP and cubic invariance is too small to be seen in comparison with the statistical errors.
Since the violation of CP and cubic invariance disappears in the large N limit, we conclude that in practice one would not be able see any effect of this symmetry breaking for even larger values of N . Hence, hereafter, we show the coefficients averaged over all µ-ν planes (S = S 1 + S 2 ) in order to study the k dependence and N dependence. We emphasize the small errors of our coefficient determinations, typically of order 10 −5 . This gives about 3 to 5 significant digits in some determinations.
We next investigate the dependence on the flux parameter k at N = 49. For that purpose we studied the k = 1 and k = 3 results at fixed value of N MD = 32. The results for the coefficients for the plaquette and 4 × 4 Wilson loop are displayed in figure 3  correspond to the analytic values. For the plaquette, the observed k-dependence is compatible with statistical errors. As expected, the differences are small but well beyond statistical errors for the 4 × 4 loop. This is also the case for the third order coefficientŴ RR 3 for which we have no analytic results. Anyhow, no dramatic k dependence is seen for the third and fourth order coefficients.
Our results averaged over all planes are displayed in table 4. The choice of values (N, k) = (16, 1), (49,2), (121, 3) are selected so as to keep k/L at ∼ 0.27 (k/L = 1/4, 2/7, 3/11 for each data set). The results at the one and two-loop level are consistent with the analytic values within two standard deviations. We also display (NPT fit) an estimate of the three-loop coefficients for N = 16 and 49 obtained by fitting nonperturbative expectation values to a third order polynomial in λ with the first coefficients fixed to the analytic values [24]. The results are roughly consistent with our NSPT determinations, which are expected to be more reliable.
Since the main interest of the TEK model is its large N limit, which coincides with Yang-Mills at infinite N and infinite volume, we give the infinite N extrapolation of our results. As mentioned previously finite N corrections are expected to grow for large values of R/L. Thus, we exclude the data of (N, k) = (16, 1) and base our analysis on (N, k) = (121, 3) and (49, 2) only and assuming a 1/N 2 dependence of the coefficients. The results are displayed in figure 4 and     would be necessary to have at least 3 sufficiently large values of N , which is left out of this exploratory work.

Analysis of the distribution
For the purpose of determining the statistical requirements involved in a precise determination of the perturbative coefficients we performed an analysis of the probability distri- bution of the corresponding NSPT estimates. As a quantitative measure we computed the cumulants up to the 5th cumulant. The N dependence of the variance (2nd cumulant) is important since the statistical error at fixed statistics is proportional to this variance. The non-zero value of higher order cumulants measures the deviation from the normal distribution. This is of practical importance since, as seen in ref. [63], the distribution for perturbative coefficients at higher order in NSPT could show strong deviations from a normal distribution, the so-called "Pepe effect". We show the histograms forŴ 11 andŴ 33 in figures 5 and 6, respectively. We have fixed N MD = 32 for better comparison. The distributions are plotted in logarithmic scale for which the normal distribution corresponds to a parabola. The dotted lines are the best fit to a normal distribution. One can see deviations in the tail for some plots for N = 16, but it disappears for larger N . Hence, the problem of "Pepe effect" is absent in the TEK model for sufficiently large N , similarly to what happens for pure SU(3) lattice gauge theory. Since the 3rd-5th cumulants have minor effect or are not detectable on the distribution, we focus on the variance only in the following analysis.
We have studied the dependence of the variance ofŴ RR on R, and N . Notice that we can also use the results for half-integer in this analysis. The R dependence is illustrated in figure 7, and is seen to go linearly with R 4 . On the other hand the dependence on j = 2 , displayed in figure 8 shows an exponential behaviour for sufficiently large j. Altogether, an ansatz of the form with c = 0.0013 (solid lines in these figures) describes our data quite well. For large N this corresponds to a dependence on the ratio of the loop size with respect to the effective size of the box R/L to the fourth power, similar to the observed subleading corrections in perturbation theory. The 1/N 2 dependence of the variance follows from the perturbative realization of factorization W 2 = W 2 + O(1/N 2 ). Similar arguments lead to the prediction that the p-th cumulant scales as O(1/N 2p−2 ). This by itself explains why the distribution tends to normal at large N . We cannot check this behaviour of the higher cumulants since their values are too small.
The exponential dependence of the variance on the order 4 coincides precisely with the characteristic growth in the number of planar diagrams [27,64,65], given by the Catalan number. In the absence of renormalons, a similar growth is also expected for the mean  valuesŴ , but it is hard to check this behaviour with our results reaching only up to order = 4. Indeed, as observed in ref. [47], extremely large order perturbation coefficients is needed to detect the expected factorial behavior for the pure SU(3) lattice gauge theory. It would be interesting to extend our results to higher order to explore this point. The computational requirements will be studied in the next subsection.
4 Estimate of the computational cost at high order and large N Here we will evaluate the computer requirements to extend our calculationŴ up to order = N trunc /2. The precision in the determination is proportional to the standard deviation of the corresponding probability distribution. Hence, if we want to keep this precision fixed our data sample should grow with the square, the variance. We previously estimated that for each sample point the computational cost goes as O(N 7/2 N 3 trunc ) for the MD part and as O(N 3 N 4 trunc ) for the reunitarization part. Using eq. (3.1), we can then give an estimate of the total computational cost as follows: This estimate can be reduced by improving the algorithm. To obtain perturbative coefficients at very higher order, the exponential scaling behaviour, 2 Ntrunc , has to be relaxed by reducing the variance using a method such as a kind of reweighting method in NSPT. Furthermore, the factor N 4 trunc comes from our implementation of the reunitarization algorithm, for which a better algorithm could be applicable. One possibility is to use the perturbative Gram-Schmidt algorithm for the reorthogonalization followed by the perturbative subtraction of the U(1) phase. The estimation of the U(1) phase still requires the perturbative trace-log computation whose cost is of O(N 3 trunc ). At least one can reduce the factor from N 4 trunc to N 3 trunc . In that case eq. (4.1) would dominate the total cost. The cost of O(N 2 trunc ) in the convolutional product of two series could be ameliorated by using the fast Fourier transformation (FFT) algorithm, which reduces the cost to O(N trunc log N trunc ). Improvement in this direction can be used in future studies.

Summary
In this paper we have developed and applied an HMD-based NSPT algorithm to compute the perturbative expansion of Wilson loops in the TEK model. At large N these coefficients coincide with those of large volume Yang-Mills theory, hence their interest. The low order coefficients are very precise and match with the values obtained by the analytic calculation of ref. [24]. We have been able to extend the perturbative series two more orders, up to O(λ 4 ).
We also studied the statistical properties of the coefficients in NSPT. We found that their distribution tends to normal at large N and there is no Pepe effect. Furthermore, we have determined the dependence of the corresponding variance on N and the order of perturbation theory. This has allowed us to estimate the computational requirements necessary to extend our results to higher orders. This would allow extending studies of the type described in ref. [27] to Yang-Mills theory in the limit of infinite number of colours. support by Priority Issue 9 to be tackled by using Post K Computer. The numerical computation for this paper were performed on the subsystem A of the ITO supercomputer system at Research Institute for Information Technology, Kyushu University, and workstations at INSAM (Institute for Nonlinear Sciences and Applied Mathematics), Hiroshima University.

A Perturbative matrix functions
In this appendix we describe the details of evaluating the perturbative coefficients of a matrix function with a perturbative series as the argument. After introducing generic algorithm for perturbative matrix functions, we explain the cases with matrix exponential and matrix logarithm explicitly.
We consider a matrix function f (F ) whose argument F is a matrix. We assume that f (F ) has the following Taylor expansion form: We also assume that the perturbative expansion of F is where g is the coupling constant and F (k) are the coefficient matrices. Note that the leading term of F is O(g). We truncate the perturbative expansion at N trunc -th order. We would like to know the coefficient [f (F )] (k) in terms of c k 's and F (k) 's. A recursion relation to evaluate the coefficient [f (F )] (k) can be obtained by expanding the Horner's method for polynomials as follows.
As F is truncated at N trunc -th order, the Taylor expansion is also truncated at N trunc -th order, where the matrix G is introduced for the truncated form. G can be evaluated via the Horner's method: G = f (F ) = c 0 + c 1 F + c 2 F 2 + · · · + c Ntrunc F Ntrunc where G j are working area. As F has the series form eq. (A.2), G j also has the series form: Here we discard higher order terms O(g k ) with k > N trunc . Consequently the perturbative coefficient [f (F )] (k) of f (F ) can be obtained as algorithm 3. The computational cost scales with N 3 trunc .

B Perturbative reunitarization
Reunitarization for U is applied perturbatively on U (k) via the perturbative expansion of the matrix logarithm: where A should be anti-Hermitian and traceless for U to be SU(N ). From the perturbative expansion of the matrix logarithm and A = ∞ k=1 g k A (k) , we have In order to apply algorithm 4 to the TEK model, we have to reunitarize the matrix V is computed.