Perturbative contributions to Wilson loops in twisted lattice boxes and reduced models

We compute the perturbative expression of Wilson loops up to order $g^4$ for SU($N$) lattice gauge theories with Wilson action on a finite box with twisted boundary conditions. Our formulas are valid for any dimension and any irreducible twist. They contain as a special case that of the 4-dimensional Twisted Eguchi-Kawai model for a symmetric twist with flux $k$. Our results allow us to analyze the finite volume corrections as a function of the flux. In particular, one can quantify the approach to volume independence at large $N$ as a function of flux $k$. The contribution of fermion fields in the adjoint representation is also analyzed.


Introduction
Within lattice gauge theory calculations, finite volume perturbative studies are interesting for various reasons, one being that numerical results are always at finite volume. From the first studies it was clear that the periodic boundary conditions (PBC) introduced complications associated to the existence of infinitely many gauge inequivalent zero-action configurations: the torons [1]. Furthermore, the toron valley is not a manifold, but rather an orbifold possessing singular faces and points. Considerable effort was put into setting up a consistent computational weak coupling expansion [2][3][4][5][6]. Other studies simply ignored the problem by expanding around a single type of minima [7]. The results should approach those obtained at infinite volume [8][9][10][11][12] 1 . 't Hooft realized that PBC are not the only possible boundary conditions for SU(N) gauge theories on the torus. He introduced the concept of twisted boundary conditions (TBC) [13,14] that was soon translated to lattice computations [15]. The new boundary conditions associate to each plane of the torus a certain flux defined modulo N , collected into an integer-valued twist tensor n µν . Already in the first studies it became clear that TBC introduced considerable simplification in perturbative calculations at finite volume [1].
The size of finite volume corrections was found to be directly connected to the magnitude of N , the number of colours of the theory. This followed from the observation made by Eguchi and Kawai [16] when studying the Schwinger-Dyson equations for Wilson loops. Their claim can be phrased as the statement that, under certain assumptions, finite volume corrections vanish in the large N limit. This leads to the large N equivalence of ordinary lattice gluodynamics with matrix models obtained by collapsing the lattice to a single point: Reduced models. If volume independence holds in the weak coupling region, this should show up in the perturbative expansion of Wilson loops. This was found to be false [17] for the original proposal of Ref. [16] (Eguchi-Kawai model). The problem arises from the attraction among the eigenvalues of the Polyakov loops induced by quantum corrections. This invalidates the Z 4 (N ) center-symmetry assumption of the equivalence proof. This could have been anticipated on the basis of the results of Ref. [1].
Identification of the source of the breakdown allowed the authors of Ref. [17] to propose a modification, called the Quenched Eguchi-Kawai (QEK) model, which could solve the problem. In this proposal the expectation values were computed taking these eigenvalues as frozen or quenched. The final results were then averaged over them. Within the perturbative regime this idea was analyzed in Ref. [18] as part of a general framework called the quenched momentum prescription. It was shown how the reduced model reproduced the perturbative expansion of the full theory. Indeed, the aforementioned eigenvalues played the role of effective momentum degrees of freedom. This particular connection between internal and space-time degrees of freedom is quite general as shown by Parisi [19].
In all the previous works periodic boundary conditions were assumed. However, two of the present authors [20] argued that Eguchi-Kawai proof holds also for TBC, which as mentioned earlier have a very different weak-coupling behaviour. This allowed them to present a reduced model, called the Twisted Eguchi-Kawai model [21] (TEK), which could achieve the large N volume independence result at all values of the coupling. With a suitable choice of the twist-tensor n µν one is guaranteed to have zero-action solutions without the zero-modes (torons) which complicate the perturbative expansion in the absence of twist. One particular simple choice of twist is the so-called symmetric twist which demands N =L 2 and has a common flux |n µν | = kL. In this case the classical vacua break the Z 4 (N ) center symmetry of the model down to Z 4 (L). This remnant symmetry is enough to guarantee the volume independence of loop equations in the large N limit.
The authors of Ref. [21] then considered the perturbative expansion of the model by expanding around any of the N 2 gauge inequivalent vacua. The Feynman rules were obtained and this iluminated the way in which the infinite volume theory is recovered from the matrix model. An important ingredient is the use of a basis of the Lie algebra of the group which has the form of a Fourier expansion. This illustrates a new and more efficient way in which space-time degrees of freedom are obtained from those in the group: the N 2 degrees of freedom of the U(N) group show up as the spatial momenta of anL 4 lattice (colour momenta). The idea can be used to achieve a volume reduction of theories with scalar or fermionic fields and can be extended to the continuum [22] 2 . A similar treatment was done for d = 2 and non-gauge theories in Ref. [23].
A bonus of this perturbative construction is that it gives a hint of what happens at finite N . The propagators are identical to lattice propagators at finite volume. Thus, finite N corrections appear in part as finite volume corrections. This is not the end of the story because the Feynman rules of the vertices adopt a peculiar form, including colourmomentum dependent phase factors. In Ref. [21], the authors showed how these phase factors cancel out in planar diagrams. It is these surviving phases that are instrumental in suppressing non-planar diagrams and reproducing the perturbative expansion of large N gauge theories. Many years later [24] the origin of these peculiar phase factors was clarified as a distinctive feature of field theories in non-commutative space-times (for a review see Ref. [25]). This led some authors [26][27][28] to propose the use of the TEK model as a regulated version of gauge theories of this type, somehow inverting the path that led to Ref. [21].
The interest of studying (lattice) gauge theories with TBC beyond the large N reduced model context was soon emphasized by several authors [29][30][31][32] and the perturbative technique extended to include space-time momenta added to the colour momenta. Since then, several perturbative calculations have been performed with different observables and contexts in mind [33][34][35][36][37][38][39].
Our present work focuses on the perturbative expansion up to order g 4 of Wilson loops for an SU(N) lattice gauge theory with Wilson action in any dimension and in a box with any irreducible orthogonal twisted boundary conditions. This means that the twist must allow the existence of discrete zero-action solutions and no zero-modes. This includes the case of the symmetric twist used in the TEK model. Our formalism is developed for any box size. In this way we bridge the gap between the infinite volume perturbative results and the L = 1 TEK model. Some preliminary results were presented in the 2016 lattice conference [40].
There are many interesting issues that our analysis aims at elucidating. These are connected to the interplay between the different parameters entering the game: the box size, the rank of the matrices N and the integer fluxes defining the twist. One of the aspects has to do with the approach to the large N limit. Volume independence would imply that in this limit the results should not depend on the lattice size. However, it is interesting to ask, as the authors of Ref. [41,42] did in the periodic boundary conditions case, what is the optimal balance of spatial and group degrees degrees of freedom that minimizes corrections. This is intimately connected to the important practical problem of estimating the corrections to volume independence for large but finite N . Depending on the results, the usefulness of reduced models as an effective simulation method to compute observables of the large N theory could be severely limited. From a more conceptual viewpoint one would like to understand the nature of these corrections. As mentioned earlier, some of these finite N corrections amount to finite volume corrections with an effective volume which depends on N ( √ N for the 4 dimensional symmetric twist). However, that is certainly not all. Some effects do depend on the phase factors at the vertices, which are a function of the twist tensor. The relevance of monitoring this dependence has been recognized recently in certain non-perturbative studies of the TEK model. At intermediate values of the coupling, several authors [43][44][45][46] reported signals that the center symmetry of the four dimensional TEK model was broken spontaneously. This is crucial, since in that case the proof of volume independence of Eguchi and Kawai fails. The problem was analyzed in Ref. [47], concluding that to avoid the problem one should scale appropriately the unique flux parameter k of the model when taking the large N limit. The validity of volume reduction under these premises has been verified in very precise measurements of Wilson loops [48]. Similar constraints are found when analyzing 2+1 dimensional theories defined on a spatial torus with twist [39,49]. In that case the problem of finding an optimal flux is related to some recent conjecture in Number Theory [50].
Obviously, all these problems do not arise in perturbation theory since centre-symmetry cannot be broken in our finite N and finite volume setting. However, the analytic calculations of perturbation theory can give hints about the origin of the possible transitions occuring when taking the volume or N to infinity. We emphasize that our computation, being of order λ 2 , already includes self-energy and vertex gluonic contributions which contain ultraviolet divergences in the continuum limit. This also relates to problems reported in the perturbative expansion of non-commutative field theories [51][52][53][54][55][56][57] having to do precisely with the self-energy of the gluon. We recall that the twisted theory can be seen as a regulated version of Yang Mills theory on the non-commutative torus. Indeed, some instabilities also arose when analyzing the model within that context [58].
The lay-out of the paper is as follows. In section 2 we set up the methodology. Our presentation is mostly self-contained and general enough to cover all the twists of the allowed type. At particular points we focus on the specific situation for symmetric twists in 2 and 4 dimensions. To facilitate the reading the Feynman rules necessary to perform the calculation are collected in two appendices. In section 3 we present our results, first to order λ ≡ g 2 N , and next to order λ 2 . These results appear as finite sums over a range of momentum values which depends on the twist tensor. In the next section (section 4) we analyze these results. In particular, we split the contributions into sets and study the difference between the computations with twist and those obtained with periodic boundary conditions ignoring the contributions of zero-modes. The focus is on the analysis of the dependence N and the box size specially when any of the two is large. For practical reasons our analysis is concentrated on the case of a symmetric box with a symmetric twist where there is only one size parameter L and one flux value k. This is specially the case in section 5 where some of the sums are evaluated numerically and the results analyzed. Comparison of the finite N corrections for the reduced model and other partially reduced options is also addressed. Several formulas that allow the analytic calculation of the leading finite volume corrections are collected in the last two appendices of the paper.
In Section 6 we try to make our analysis more complete by analyzing a few extensions. In particular we consider the distinction between U(N) and SU(N) results and the additional contributions to the Wilson loop coming from the quarks in the adjoint representation of the group. The latter can be included in a twisted setting in a rather straightforward way, while quarks in the fundamental need the addition of replicas (flavours). The explicit formulas for fermions have been obtained for Wilson fermions at r = 1 and critical value of the hopping. In that same section we also compare our results with high statistics measurements of the loops with standard Monte Carlo techniques. Our goal is to be able to test extremely tiny effects such as the breakdown of cubic invariance by the twist as well as the non-zero value of the imaginary part. The data match perfectly with the expectations. Furthermore, this analysis also allows for an estimate of the coefficient of order λ 3 and the determination of the range of couplings for which the truncated perturbative expansion is a good approximation.
The paper closes with a conclusions section in which we sum up the main results following from our calculation.

Models and methodology
In this section we will describe the type of models that we will be considering as well as the tools necessary for the calculation of the coefficients.

The action with twist
We will be considering d-dimensional SU(N) lattice gauge theory with Wilson action on an hypercubic box of size L 0 × · · · × L d−1 . The L µ can be taken as the components of a d-dimensional vector L. The product of the sizes gives the total lattice volume labelled V . Different lengths in different directions break the hypercubic invariance. Hence, for simplicity we will often specify results for a symmetric box L µ = L. The action depends on a single coupling b = β/(2N 2 ). We will focus upon the behaviour of the expectation values of R × T rectangular Wilson loops: where U (R, T ) is the ordered product of all links around the perimeter of the R × T rectangle. Our goal is to study the behaviour of these observables for large values of b. In this limit the result is calculable using perturbative/weak-coupling techniques. As mentioned in the introduction, this problem has been addressed earlier by several authors [7][8][9][10][11][12]. The main difference of our work with others is that we will consider arbitrary orthogonal irreducible twisted boundary conditions on the lattice [13,59]. In particular this will include symmetric twisted boundary conditions in four dimensions [21]. We do not intend to review the formalism to implement twisted boundary conditions on the lattice [15,60]. We will just remind the readers that after a change of variables on the links one reaches an action of the form where the link matrices are periodic U µ (n) = U µ (n + L νν ), and the plaquette factors Z µν (n) = Z * νµ (n) are elements of the center. Not all values of Z µν (n) amount to twisted boundary conditions. First of all, it is necessary that the product of all Z factors over the faces of every cube (taken with orientation) is equal one. The non-trivial twist follows by multiplying all the Z factors in each µ − ν plane, to give an overall center-elementẐ µν : Because of the condition on cubes, theẐ µν do not depend on the position of the plane but only on its orientation. Being elements of the center one can writeẐ µν = exp{2πin µν /N }, where n µν is an antisymmetric tensor of integers defined modulo N . It is this twist tensor that specifies the twist. Redefining the link variables by multiplication with an element of the centre, one can change the value of the individual plaquette factors Z µν (n), but the twist tensor remains unaffected. This allows one to set all the plaquette factors to 1, except for a single twisted plaquette in each µ − ν plane. Conventionally, one can choose that plaquette to be one at the corner (n µ = L µ − 1; n ν = L ν − 1). The change of variables that led to the action Eq. (2.2), also transforms the Wilson loop expectation variables. These are modified as follows where Z(R, T ) is the product of the Z µν (n) factors for all plaquettes which fill up the rectangle. Notice that the result will also depend on the directions defining the plane in which the rectangle is sitting.
In the next subsections we will derive the perturbative expansion of these quantities up to order 1/b 2 . We will specify the meaning of orthogonal irreducible twists and provide some examples in various space-time dimensions.

Classical minima of the action
As b −→ ∞ the functional integral is dominated by the configurations that minimize the action. We will restrict ourselves to twist tensors for which the corresponding minimum action vanishes. These are called orthogonal twists 3 . The corresponding zero-action configuration will be named as follows U µ (n) → Γ µ (n). The zero-action condition implies that the SU(N) matrices Γ µ (n) satisfy Γ µ (n)Γ ν (n +μ) = Z νµ (n)Γ ν (n)Γ µ (n +ν) (2.5) Obviously, the solution is not unique, since any gauge transformation of this one gives also a zero-action solution: where Ω(n) are arbitrary SU(N) matrices periodic on the torus. New solutions can also be obtained by the replacement where z µ is an element of the center. In some cases these solutions are gauge inequivalent to the previous ones. To study this point one must analyze the remaining gauge invariant observables of this zero-action configurations: the Polyakov loops. For every lattice path γ with origin in one lattice point, which we label 0, and endpoint n, one can construct the path-ordered product for this zero-action configuration which we will label Γ(γ). The closed lattice paths can be classified into subsets according to its winding numbers around each of the torus directions. The condition Eq. (2.5) implies that the contractible loops (which are associated to vanishing winding) are given by elements of the center. Similarly, those paths having the same number of windings have path-ordered exponentials which are unique up to multiplication by an element of the center. Now we can choose one representative path having winding 1 in the direction µ and zero in the remaining torus directions, we will call its associated path-ordered exponential Γ µ . From Eq. (2.5) we deduce that these SU(N) matrices must satisfy whereẐ νµ are the elements of Z N introduced earlier and characterizing the twist. Notice that the previous equation does not depend on the choice of representative paths. Indeed, it is the existence of solutions to Eq. (2.8) what guarantees the existence of zero-action solutions and the definition of orthogonal twists. For a more thorough discussion of the conditions on the twist tensor n µν we refer the reader to Ref. [59]. Furthermore, we will restrict ourselves to what we call irreducible twists. These are defined by the equivalent of Schur lemma, stated by saying that the only matrices which commute with all Γ µ are the multiples of the identity. If we restrict ourselves to SU(N) matrices, the set of gaugeinequivalent solutions becomes discrete. From a practical viewpoint, the irreducibility condition eliminates the presence of zero-modes which complicate the perturbative analysis considerably. Hence, the Γ µ matrices are uniquely defined up to a unitary similarity transformation, which is just a global gauge transformation, and a multiplication by an element of the centre. The second operation defines the centre symmetry group. However, not all these transformations produce gauge-inequivalent solutions. This only occurs when the eigenvalues of the Γ µ matrices, which are gauge invariant quantities, change [59].
The irreducibility condition implies that the algebra generated by multiplication of the Γ µ matrices has complex dimension N 2 . It also implies that Γ N µ must be a multiple of the identity for all µ.

Derivation of the perturbative expansion
To proceed with the perturbative expansion one has to expand the link matrices around the zero-action solutions as follows This is just an expansion around a background field and, as is well-known (see for example [61,62]), the plaquette becomes where we have introduced the forward lattice derivative ∇ + µ . In our case it is actually a covariant derivative with respect to the background field given by the zero-action solution. The expression of G µν (n), obtained by applying the Baker-Campbell-Haussdorf formula, is similar to the one obtained for periodic boundary conditions with the primes modifying the translated vector potential. For example, the leading term is The lattice vector potentials, for each link direction ρ, are V traceless hermitian N × N matrices (V is lattice volume). This is a V (N 2 − 1)-dimensional real vector space, and we can take as its basis the simultaneous eigenstates of the ∇ + µ operators (notice that they commute). We call these basis vectors χ(n; q) (which are not necessarily hermitian) and they satisfy ∇ + µ χ(n; q) = (e iqµ − 1)χ(n; q) (2.13) The form of the eigenvalues comes from the definition of the operator ∇ + µ . Spelling out the condition one must have Γ µ (n)χ(n +μ; q)Γ † µ (n) = e iqµ χ(n; q) (2.14) To solve this equation we choose one reference point on the lattice, which without loss of generality we fix as n = 0. For any other point we choose a non-winding forward moving path γ(n) joining the origin with that point. Then we have χ(n; q) = e iqn Γ † (γ(n))χ(0; q)Γ(γ(n)) (2.15) Notice that the solution does not depend on the choice of path γ(n), because the corresponding matrices differ by multiplication by an element of the center. The final requirement is that the eigenvectors satisfy the required periodic boundary conditions. Defininĝ Γ(q) ≡ χ(0; q), this condition implies that for any direction µ we must have This is a well-studied matrix equation. The condition of irreducibility implies that there are (N 2 − 1) traceless linearly independent solutions. From irreducibility one can conclude that L µ q µ must be an integer multiple of 2π/N . Hence, we can write L µ q µ = 2πmµ N , where the integers m µ are defined modulo N . If we remove the condition of vanishing trace there is an additional solution given by a multiple of the identity matrix and having m µ = 0. Now, coming back to the original eigenvalue equation Eq. (2.14), and realizing that the q µ are defined modulo 2π, we conclude that we have a total of V (N 2 − 1) different eigenvalues, each characterized by a different d-dimensional vector q. These momentum vectors have the form where the m µ are the integers introduced earlier, which enter in the first term which we call colour momentum. The second term has the standard form of momenta in a periodic lattice and is thus labelled spatial momentum. It is convenient to include also the m µ = 0 solution because then the set of momenta has the structure of a finite abelian group, which we will call Λ mom . It is a subgroup of the group Furthermore, the set of spatial momenta Λ L is a subgroup of Λ mom having V elements. Colour momenta are more rigorously identified with elements of the quotient group Λ mom /Λ L , having N 2 elements. Let us now focus on the eigenvector χ(n; q). The eigenvalue equation only fixes these matrices up to multiplication by a constant. Part of the arbitrarity can be fixed by a normalization condition. We will impose 1 V n Tr(χ † (n; p)χ(n; q)) = 1 2 δ p,q (2.18) which fixesΓ(q) to be a unitary matrix divided by √ 2N . This leaves a phase arbitrarity which can be further reduced by imposing additional conditions on the unitary matrix. For example, one can impose that it belongs to SU(N). Alternatively, we can impose that ( √ 2NΓ(q)) N = I. Any of the two conditions, which might be incompatible with each other as we will see later, reduces the arbitrarity to multiplication by an element of the center Z N . A choice of this element for every q fixes the assignment q −→Γ(q) (2.19) This provides a group homomorphism from Λ mom to SU (N )/Z N . For the SU(N) normalization condition, this impliesΓ where Φ(q, p) is an integer multiple of 2π/N , which depends on the choice of phases. We can restrictΓ(0) by demanding Φ(0, q) = Φ(q, 0) = 0. If we adopt the ( √ 2NΓ(q)) N = I condition, Eq. (2.20) also holds, but then e iΦ(q,p) could be an element of Z 2N .
Having solved the eigenvalue and eigenvector problem, we realize that given that our solutions are a collection of linearly independent matrix fields, we can actually decompose our vector potentials as follows which generalizes the Fourier decomposition. Two comments are necessary at this point. The first affects the condition that the vector potentials are hermitian matrices. This imposes a constraint on the Fourier coefficientsÂ µ (q) as follows: The second is the requirement of tracelessness, specific of SU(N). This implies thatÂ(0) = 0 for q ∈ Λ L . Hence, the sum extends over the set difference Λ mom \ Λ L . For simplicity this restriction will be noted by the prime affecting the summation symbol.
Combining the previous expression we can define which play the role of the d and f symbols of the SU(N) Lie algebra in our basis. By definition D is completely symmetric and F completely antisymmetric under the exchange of their arguments. What is the connection between the choice of twist tensor and the value of the lattice of momenta Λ mom ? What is the explicit form of the matricesΓ(q) and of the F and D functions? This can be analysed as follows. As mentioned previously the matrices Γ µ generate an algebra by multiplication. By irreducibility, this algebra has dimension N 2 . Hence, the matricesΓ(q) must necessarily have the following form where α(q) are integer multiples of π/N and s µ (q) are integers defined modulo 2N . Using the commutation relations of the Γ µ one can find the relation between the integers s µ (q) and q, given by The previous equation defines a homomorphism N from the group (Z/2N Z) d to Λ mom /Λ L . This cannot be an isomorphism except in two dimensions since the number of elements in Λ mom /Λ L is N 2 , which is smaller that (2N ) d . Indeed, using the isomorphism theorem we conclude that This allows the computation of Λ mom given the twist tensor. On the other hand the inverse q −→ s(q) is not uniquely defined and is a matter of convention. This convention dependence is also present in the choice of elements α(q). Indeed, any choice of inverse can always be compensated by appropriately choosing the α. The convention dependence extends to the value of Φ(q, p) and the F and D symbols. A convenient choice is to impose the condition Φ(p, −p) = 0. This equation makes the hermiticity condition Eq. (2.22) look just like in the ordinary Fourier decomposition of a real field. It should be noted though, that for even values of N the condition might conflict with the SU(N) normalization condition. In combination with the alternative condition ( √ 2NΓ(q)) N = I, it fixes the value ofΓ(q) up to a sign. We stress, nonetheless, that the convention adopted for the definition ofΓ(q) affects only the corresponding definition of the Fourier coefficientsÂ µ (q) and has no influence in the results of Wilson loops or other observables.
Furthermore, it is important to realize that the antisymmetric combination of the phases (sum over repeated indices implied) is convention independent. The previous equation is an equality among angles, and hence is defined modulo 2π. The antisymmetric matrixñ µν is defined by the relation n µαñαβ n βν = n µν mod N . Its matrix elements are not necessarily integers, but the inversion formula (see below) should be well-defined. Its existence can be deduced by transforming n µν to its canonical form (see Ref. [59]). Although the matrix is not unique, its arbitrarity does not affect Eq. (2.27). Its non-uniqueness however shows up when using the matrix to define the inverse map q −→ s(q) as follows The condition that the left-hand side are integers provides the restriction on the elements n µν , mentioned above. To summarize, we can say that up to now the presentation has been completely general for the case of irreducible twists 4 . The most important ingredients are the presence of a lattice of momenta Λ mom and the convention dependent value of the D and F symbols. In the next subsection we will apply our formalism to the most useful cases in two, three and four space-time dimensions, and provide explicit formulas for the different ingredients in terms of the twist tensor. The four-dimensional case is the most interesting and will be used for most of the numerical analysis that will follow later.

Particular cases of twists in 2 to 4 dimensions
The two dimensional case is particularly simple since twist tensors are of the form n µν = k µν , where 01 = − 10 = 1. The condition of irreducibility amounts to constraining the integer k to be coprime with N . The lattice of momenta Λ mom is simply given by all momenta having the form q µ = 2πmµ LµN , where m µ are integers modulo L µ N . Notice that this is equivalent to the standard lattice momenta in a box of size (N L 0 ) × (N L 1 ), with an effective volume of V eff = N 2 V . The Γ µ matrices can be written in terms of 't Hooft matrices Q and P satisfying P Q = zQP (2.29) where z = exp{2πi/N }. The matrices are given by Q = diag(1, z, z 2 , . . . , z N −1 ) and P ij = δ j i+1 . Thus a possible choice of matrices satisfying the algebra would be Γ 0 = Q and Γ 1 = P k . Notice, however, that for even N the matrices have determinant −1. Thus, if we impose that Γ 0 belongs to SU(N) one should rather take Γ 0 = ±iQ and Γ 1 = (±iP ) k , but paying the price that now Γ N 0 = −1. This is the conflict of normalization conditions that we were mentioning earlier. The same normalization conflict translates to the choice ofΓ(q). We might obviously writê Herek is an integer defined through the relation: We only need to fix the function α(q). One way to fix it is to demand that the hermiticity condition Eq, (2.22) adopts the same form as for ordinary Fourier expansion, namely setting Φ(q, −q) = 0. This leads to where S(m 0 , m 1 ) takes the value 0 or 1, giving the two possible values of the square root. This second term is necessary because it compensates for the fact that the first term is not always invariant under the shift m µ −→ m µ + N L µ . We might set it to zero if we restrict m µ to lie in a particular range. If N is odd one can always choosek to be even (ifk was odd, replace it byk + N ) and one can directly set S(m 0 , m 1 ) = 0. In that case one finds For the even N case, the formula is still valid if we set the integers m µ to lie in a particular interval.
For the three dimensional case the twist tensor can be written in terms of the completely antisymmetric symbol with three indices as follows: n µν = µνρ r ρ , where r is a vector of integers modulo N . The irreducibility conditions amounts to the fact that the greatest common divisor or r α and N is 1. The integers s µ (q) must satisfy where we have used the standard three-dimensional notation for vector products. Hence, the space of momenta Λ mom can be identified with those that correspond to m· r = 0 mod N . The irreducibility condition now guarantees that there exist a vector of integers v such that r · v = −1 mod N . This allows to define a possible inversion as s(q) = v × m. The rest follows similarly to the two dimensional case withñ µν = µνρ v ρ . In four dimensions, the orthogonality condition requires a twist satisfying κ(n µν ) ≡ µνρσ n µν n ρσ /8 = 0 (mod N ). Irreducibility is granted provided the greatest common divisor of N , n µν , and κ(n µν )/N is equal to 1. The case in which the twist is in one plane or in a three dimensional section proceeds identically to the previous cases with q c living in a 2 or 3 dimensional subspace as before. Hence, we will now focus in the case of the symmetric twists where N is the square of an integer N =L 2 and The twist is irreducible if k andL are coprime integers. The lattice of momenta Λ mom is given by q µ = 2πm µ /(L µL ), with m µ integers defined modulo L µL . This leads, like in the two dimensional case, to an effective lattice volume V eff = N 2 V . The integers s(q) are given by Herek is an integer defined through the relation: and˜ µν is an antisymmetric tensor satisfying: Notice that this defines theñ µν matrix to be given byk˜ µν /L. The function Φ(p, q) becomes: As in the two dimensional case, imposing hermiticity by setting Φ(q, −q) = 0 leads to: A particular choice satisfying this condition for a momentum q µ = 2πm µ /(L µL ) is: leading to D and F functions defined in terms of θ(p, q) as in Eq. (2.35), with: where we have introduced the antisymmetric tensor θ µν defined as: with the angleθ ≡ 2πk/L.

The gauge fixed action at order λ
We will be using the standard covariant gauge fixing term with gauge parameter ξ. Its contribution to the action is where ∇ − µ is now minus the adjoint of ∇ + µ . This is a typical background field gauge. To order g 2 the ghost action corresponding to this gauge fixing is given by: The ghost fields c-c have a similar colour-space Fourier decomposition as the gauge fields. There is also an additional a contribution to the action coming from the expression of the Haar measure on the group in terms of integration over the Fourier coefficientsÂ(q). To order λ it is given by To derive this expression we parameterize a generic group matrix element as in terms of the basis of the SU (N ) algebra given byΓ(q c ) with q c the colour momentum taking N 2 − 1 values. The prime over the sum indicates that zero momentum is excluded. The volume element of the group in terms of these variables is with the metric G defined as: Inserting the expression for the U matrices leads to and from here one computes the O(λ) contribution to the action given by (2.53) In obtaining this result we have used the hermiticity relation on the coefficients w(q c ) and the equality which is the expression of the quadratic Casimir in the adjoint representation in our basis (see Appendix D). Summarizing, we have obtained the gauge fixed partition function to order g 2 given by This action can be expanded in powers of g to derive the Feynman rules of the theory. For example, in Feynman gauge the propagator of the gauge field reads: 56) and the ghost propagator is Notice that if we adopt the hermiticity condition on theΓ(q), giving Φ(q, −q) = 0, the expressions simplify. In what follows we will adopt this convention. Because of the problems associated with this convention and explained in subsection 2.4, the momenta are now defined in a range and not modulo 2π. Correspondingly the momentum conservation delta functions are now strict and not modulo 2π. In any case these difficulties just affect intermediate expressions and not to the final results.

Expansion of the Wilson loops
An essential ingredient in the calculation is the expansion of the Wilson loop in powers of the vector potentials A µ (n). This expansion for the particular case of the plaquette is also necessary to derive the non-quadratic terms in the Wilson action giving the vertices of the theory.
We recall the definition of our observable given in Eq. (2.4). To process the right-hand side we replace the links by the expression Eq. (2.9). In simplified notation this gives rise to where we used the label l to represent a link instead of the conventional (n, µ) combination. The product on the right hand side is the ordered product of the exponentials around the rectangle R. Finally A l is given by where Γ l is the product of the Γ µ (n) factors from one reference point in the square to the origin of the link l. The dependence on the choice of reference point drops out when taking the trace. Notice that in terms of the A the Z(R, T ) factor has disappeared from the right-hand side.
One can now use the Baker-Campbell-Haussdorf formula to rewrite this as: with G a hermitian matrix which can be expanded in powers of g as follows: where: In the previous formula the ordering of the links is done along the perimeter of the rectangle following the plaquette orientation. Now we can express the trace in terms of G as follows: To perform the calculation we need to substitute the expression for G and use the Fourier decomposition written in the following simplified form where n are the coordinates of the lowest vertex of the rectangle. Notice that if the link l has origin n and direction µ, the coefficients are given by We will also use the following notation After averaging over n and expressing the traces in terms of the group constants F and D (using for simplicity the hermiticity condition Φ(q, −q) = 0), we arrive at: where: The previous formulas express the expectation values of Wilson loops in terms of the n-point Green functions of the vector potentials. The latter can be computed as a power series in g using the Feynman rules of the theory, given in App. A.

Results of the perturbative expansion of Wilson loops
In the present section we use the machinery developped in the previous section to compute the perturbative expansion of the expectation values of rectangular Wilson loops. In particular, we consider the coefficients of the expansion up to order λ 2 = 1/b 2 as follows: Alternatively, we might consider the expansion of the logarithm instead The two sets of coefficients are related as follows To obtain these coefficients we start by the expressions given in the previous section and expand the U (n) and V (n) terms in powers of g: Using this terminology, similar to that followed in Ref. [9][10][11], we arrive at the following expression for the coefficients of the logarithm of the Wilson loop at O(g 4 ): Notice that coefficient of order λ 2 requires the calculation at one-loop of the two point function U 2 . In the following subsections we will spell out the calculation of these coefficients.

The Wilson loop at O(λ)
Combining the previous results we obtain the expression of the first coefficient as follows: To compute this expression we must write down explicitly the expression ofĀ(q) = l∈R A l (q) in terms of the Fourier coefficientsÂ ρ (q). To simplify notation we will specify that the rectangle is sitting in the µ − ν plane, with R and T being the length of the edges in the µ and ν directions respectively. We can then separateĀ(q) as a sum of the contributions of its four edges. Noting these contributions as A (i) with i = 1, 3 (µ direction) and i = 2, 4 (ν direction), we get: where we have introduced the symbols Q µ (q) and Q ν (q) given implicitly in terms of finite geometric sums. Performing these sums explicitly we have and q µ the lattice momentum introduced in Eq. (2.58). This expression is singular for q µ = 0 in which case the result if Q µ (0) = R. Replacing µ by ν and R by T , we get the remaining symbols.
With this notation we finally get which can be rewritten in a more symmetric fashion as with: Using the expression for the propagator, gives the final result: The result agrees with the tree level result for the standard Wilson action on an infinite lattice derived in [9] if one replaces appropriately the momentum sums by integrals.
For the particular case of the plaquette (R = T = 1) the result simplifies and we get The last equality is true for the average of the plaquette over all µ−ν planes in d space-time dimensions. It coincides with the plaquette in each plane if there is symmetry among all directions. Otherwise the plaquette expectation value at this order depends on the plane.

The Wilson loop at O(λ 2 )
To compute the coefficient of the logarithm of the Wilson loop expectation value to the next orderW (R×T ) 2 (N, L, n µν ), we need to evaluate U (n) , n = 1, · · · , 4, and V (n) , n = 3, 4 in Eqs. (2.71)-(2.75) and substitute them in expression (3.6). The computation of the different terms can be done using the Feynman rules given in App. A.
In the following paragraphs we list the expression of the U (n) a and V (n) a terms entering in Eq. (3.6). As we did at leading order, we use the label µ to indicate the direction of the loop having length R and ν that having length T . We also use the simplifying symbols given below: We arrive at: The expression for the vacuum polarization Π αβ can be found in App. A.
The corresponding expressions for the plaquette simplify considerably: Notice that all the previous expressions are valid for arbitrary space-time dimension and for an arbitrary irreducible twist. All the dependence on the twist is contained in the F and D factors and in the ranges of the momentum sums. Notice also that in some of these sums we have dropped the prime affecting the summation symbol. As explained below, this is because the F factor vanishes for the excluded momenta in the primed summation.

Analysis of the results
In the previous section we have presented the result of the calculation expressed as single and double sums over discrete momenta. To help in understanding the implications it is interesting to analyze the N and L dependence and to understand the connection with the case of periodic boundary conditions. For the case of the real part, we can naturally separate out two types of contributions to the coefficients: those proportional to the structure constant square F 2 and those that are not. They will be called non-abelian and abelian respectively. The imaginary parts are always proportional to F (indeed they are also proportional to the anticommutator D) so they can be classified also as non-abelian.
We recall that the momentum sums range over a finite lattice labelled Λ mom \ Λ L , where Λ mom is a finite abelian group and Λ L the subgroup of spatial momenta. The zero momentum q = 0 (neutral element) is contained in both sets. It is convenient to exclude it from both and use a prime symbol to label the resulting set: Λ mom = Λ mom − {0}. Notice that this restriction does not affect the set difference: Λ mom \ Λ L = Λ mom \ Λ L . The removal of the zero-momentum from the sum eliminates the apparent ill-definition of the expressions. An interesting observation for the analysis that follows is that F vanishes when any of its arguments belongs to Λ L . Thus, in all contributions of non-abelian type, we can drop the prime in the sum and extend the sums to Λ mom . Now we are in position to discuss the relation between our results and those obtained for periodic boundary conditions, more precisely with the finite volume periodic results obtained by Heller and Karsch [7] by neglecting the contribution of zero-modes and explicitly excluding zero momentum in the sums. According to our previous considerations, the integrands of the different contributions to the real part of the Wilson loops for periodic and twisted boundary conditions are identical. The main difference is that for the periodic case, the momentum sums are now over Λ L , and the colour sums are performed independently. Hence, it is possible to transform our formulas into those of Ref [7] by the following substitutions: for the non-abelian terms proportional to F 2 , and for the abelian ones. In the previous formulas we have added a third column corresponding to the infinite volume limit. It reproduces the results of Weisz, Wetzel and Wohlert [11] for the four dimensional case. For simplicity let us now focus specifically in the case of a symmetric box of length L µ = L and upon symmetric twists for which all directions appear on symmetrically on Λ mom (in more detail for the d=2 and 4 cases). At leading order, the coefficient corresponding to periodic boundary conditions (PBC) is given bỹ which spells out its dependence on N and L. The function F 1 (L) is given by a single sum over momenta in the set Λ L . It vanishes for the one-point box F 1 (1) = 0. The case of two dimensions for R, T ≤ L is particularly simple since the sums can be evaluated exactly and they give For general dimension d and at large L, F 1 (L) behaves as (see Appendix C) For the particular case of the plaquette, higher order corrections vanish and the d-dimensional result is given by F 1 (L) = (1 − 1/L d )/2d. Using the formulas given earlier the result for twisted boundary conditions and symmetric twist can be expressed in terms of the same function as follows where L eff = LL is the effective size parameter, withL = √ N in 4 dimensions and equal to N in two. It is interesting to observe that the result for the TEK model (L = 1) is just F 1 (L), the large N result on a box of sizeL d . In particular, for large L the d = 4 twisted result approaches infinite volume limit value with corrections that go like (N − 1)/(N 3 L 6 ).
A similar analysis can be done at the next order. For that purpose we have to separate the different contributions into those that we called abelian and non-abelian. Within the former there are two different N dependencies corresponding to the measure term Eq. (B.1) and the abelian contribution from the tadpole Eq. (B.7) respectively (corresponding to Π mes and Π W 2 in App. B). Finally, the (PBC) result (zero-mode contribution excluded) can be expressed as followsW in terms of two functions of L. The first function F 2 (L) can be split as follows where F N A (L) includes the non-abelian part and F mes (L) the contribution from the measure. The other function F W (L) comes from another abelian contribution to the vacuum polarization and, using the symmetry of all directions, can be expressed in terms of F 1 (L) as follows: Our functions can be connected to those given by Heller and Karsch [7] for 4 dimensions by the following relations:W Now we will analyse the expression for the symmetric twist case. The different contributions can be split as follows: where the first two terms contain the same two functions that enter the periodic formula and the last two are specific of the twisted case.
where A(p, q) is a smooth function of its arguments, whose explicit form can be read out from our formulas. The important part of the previous chain of equations is the substitution N F 2 = 1 − cos, which leads to its decomposition into two functions. The first F N A was already present in the periodic case. The new function F N P (N, L, k) is the only one in which the N , L, and k dependence are mixed up. The analysis of all diagrams that was carried out in Ref. [21] implies that this term contains the contribution of non-planar diagrams, and this explains the name given to it. Notice that for volume reduction to hold in perturbation theory F N P (N, L, k) should go to zero in the large N limit. In appendix D we analyze the behaviour of this type of sums as a function of N and L. In particular we show that when L goes to infinity F N P (N, L, k) tends towards − 1 which then goes to zero when L goes to infinity. The added piece is substracted out and combined with the −1/N 2 F mes term to produce our final expression Eq. (4.15). The usefulness of making this arrangement goes beyond this simplicity. Indeed, the new function F 2T (N, L, k) contains all the twist dependence and goes to zero in the infinite volume limit and in the large N limit. The imaginary part of the coefficient has been collected into the function G 2T (N, L, k) which has no periodic counterpart. Its presence is due to a violation of CP symmetry induced by the twist vector. Obviously, it vanishes for k = 0 as well as in the infinite volume limit. Volume independence implies that it should also vanish in the large N limit.
Our goal is to use the decomposition Eq. (4.15), to analyse the N and L dependence of the coefficients. One particular case is that of the TEK model (L = 1). The formula simplifies since F i (1) vanishes. The resulting expression is quite appealing It means that the TEK coefficient is equal to the periodic one at large N computed at an effective lattice size ofL d plus an additional complex contribution coming from non-planar diagrams. The two terms correspond nicely with the two main effects of the Feynman rules of the TEK model: a propagator equivalent to that ofL d lattice and a modified vertex which affects the non-planar diagrams only. In the general case, as mentioned earlier, if the non-planar term F N P (N, 1, k) + iG 2T (N, 1, k) goes to zero in the large N limit, one recovers the volume reduction result: The large N limit of the twisted theory coincides with the infinite volume large N result. On the contrary, it is very clear that reduction does not work for periodic boundary conditions. In the large N limit the coefficient is given byW which is still size dependent. We emphasize nevertheless that what we call PBC is not the correct weak coupling expression for periodic boundary conditions. Our calculation does not take into account zero-modes. This is known to leading order, but not to the 1/b 2 order that we are calculating.
On the other extreme we can consider the behaviour of the perturbative coefficients for large values of L. The infinite volume coefficients coincide for periodic and twisted boundary conditions provided F 2T + iG 2T vanishes at large L (see appendix D). This result was to be expected, since at infinite volume boundary conditions should not matter.
A complete analytic study of the approach to large L is hard to do due to the complicated structure of the coefficients A(p, q). Functions involving single sums can be easily analysed though along the guidelines given in appendix C. Apart from the behaviour of F 1 (L) given earlier, we also study the L-dependence of F mes (L) for large L. Using the formulas developped in appendix C we show that in four dimensions the leading correction is given by: In two dimensions F mes (L) diverges logarithmically with L as: The complicated structure of F N A (L) involving double sums has prevented us from obtaining its expansion in inverse powers of L. Numerically it seems that, to a high precision, the leading 1/L 2 correction is equal and opposite to that of F mes (L). This is presumably associated to the vanishing of the vacuum polarization at zero momentum which occurs through a similar cancellation. This is another reason for expressing the results in terms of F 2 rather than separating out its two components. In summary, in the four dimensional case the function F 2 can be fitted at large L to a functional form A similar fit for the case of the plaquette was also advocated by other authors earlier [63]. The explicit R 2 T 2 dependence is not intended to be exact, i.e. γ 2 and γ 2 can also depend on R and T . However, our best fit values for square loops to be presented later on in table 2, give values which are of similar size.
LOOP  Table 1: Values at infinite volume of the functions F 1 and F 2 defined in the text. The next three columns are combinations of these numbers giving the second order coefficients W 2 andŴ 2 at large N , as well as the parameter K of Ref. [12] Notice that in the perturbative coefficients for twisted boundary conditions the functions F i (L) appear in the combination F i (L eff )−F i (L)/N 2 . This cancels the 1/V subleading terms but not the ones containing a logarithm.
Computing the large L behaviour of the functions F 2T (N, L, k) and G 2T (N, L, k) is even more difficult, beyond the fact that they should vanish both in the large N limit and in the large L limit. The situation is analyzed in appendix D. A numerical study will be presented in the next section for the four-dimensional case with symmetric twist.

Numerical evaluation of the coefficients in four dimensions
In parallel with the analysis performed in the previous section, it is interesting to study the numerical values of the perturbative coefficients for some values of the parameters (N , L and twist k). To obtain the numbers one has to perform the momentum sums. These are finite sums (for finite L and N ) that are encoded in the functions given in the previous section. The leading order coefficient depends on the function F 1 (L), expressible as single four-momentum sum. To the next order we need the functions F 2 , F 2T and G 2T . The first one is the sum of two terms: F mes (L) which is also given by a single four-momentum sum, and F N A (L) given by a double four-momentum sum instead. This involves L 8 sums, limiting the maximum value of L that can be achieved. Some numerical values were obtained in Ref. [7]. The functions F 2T and G 2T which are specific of the twisted case, share the same difficulty plus the additional one of depending on several variables: N , L and k.
Let us start by presenting our results for F i (L). As mentioned earlier F 1 (L) can be computed for large values of L (∼ 100). Since the leading coefficients in inverse powers of 1/L 2 are known analytically, one can extrapolate the results to infinite volume with very high precision. The values of F 1 (∞) are given in table 1 for square Wilson loops up to 4 × 4. In the case of F 2 (L) we have been able to compute it up to L = 34. Going slightly beyond this point is feasible but unnecessary. As mentioned earlier the behaviour for large L is well fitted by Eq. (4.23). In the case of larger loops, stable fits require the inclusion of 1/L 6 and 1/L 8 terms. The infinite volume value is given in table 1. Errors are obtained from the variation of the parameters with the fitting range. As a example of the quality  of the fit we display L 4 (F 2 (L) − F 2 (∞)) for a 2 × 2 loop in Fig. 1. Notice that while the infinite volume coefficients F i (∞) grow moderately in size with R, the leading correction coefficient goes rather like R 4 . This is apparent from the similar magnitude of γ 2 and γ 2 for all R as seen in table 2.
Using F 1 (∞) and F 2 (∞) we can compute the infinite volume perturbative coefficients at any N . The values at N = ∞ of the second coefficient in the expansion of the Wilson loop expectation value and its logarithm (Ŵ 2 andW 2 respectively) are also given in table 1. We also add the coefficient K(R, R) used in Ref. [11]. Our calculations are consistent with the precise results of Ref. [12] for the plaquette and improve by many significant digits the published results for larger loops [7,11,64]. Now let us proceed to study the new functions F 2T (N, L, k) and G 2T (N, L, k), which  In a recent work [39], the present authors advocated that physical results for SU(N) gauge theories on twisted boxes depend on these variables only through the combinations L eff = LL andθ = 2πk/L. This applies rather well to the 2+1 dimensional case both in perturbation theory and non-perturbatively [39,49,65] and to the non-perturbative calculation of the twisted gradient flow running coupling in SU (∞) [66].
The previous observation suggests that we display the functions multiplied by the effective volume V eff = L 4 N 2 versusk/L. All functions have a similar behaviour so that we will focus on F 2T for the plaquette for a plane in S 1 . This is given in Fig. 2. Different symbols describe the different values of the independent argumentsk and L. The plot contains a lot of information that we will now spell out. First of all, the data does not show any growth with rising L eff at fixed values ofk/L. This is very important since it validates the two main expectations of our previous discussion: that the function F 2T goes   02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Figure 4: Behaviour of F 2T for small values ofk/L to zero when either L or N go to infinity. Furthermore, it tells us that when the limit is taken at fixedθ the approach to zero goes roughly as 1/V eff . We cannot exclude logarithmic or other mild dependencies, but this would hardly change the conclusion. The result can be easily confirmed by studying the L dependence of the values at fixedk and N . Our data at N = 4, 9, 16, 25, 49 cover a sufficiently large number of L values to get a good fit to a 1/L 4 dependence (see Fig. 3). Concerning the N dependence the test is complicated by the fact that when we changê L we are also changingk/L. However, as we slightly change the value ofk/L the value changes only by factors of 2 or so. It is unclear at this stage whether as N gets larger one approaches a smooth oscillatory function or not. In any case, these changes are small compared to the large changes in values of V eff . Indeed, the value of F 2T itself at neighbouring points sometimes changes by three orders of magnitude. As an example, let us discuss the results for the rangek/L ∈ [0.27, 0.3]. We have 13 different values ofk , L andL which give data in this region. The values of F 2T themselves change considerably within this set. The result for L = 1,L = 7,k = 2 is 2.28 10 −6 , which multiplied by the effective volume gives 0.0055. On the other extreme we have values as low as 7.42 10 −9 , 4.05 10 −9 , 6.71 10 −9 for (L, L,k) values of (29,1,8), (17,2,5) and (7,4,2), which multiplied by the effective volume give 0.0053, 0.0054 and 0.0041 respectively. Similar results (within a factor of 2) are obtained for the remaining data points. We believe this is enough to put our main conclusion on robust grounds.
A different perspective is obtained if instead of fixingθ we fixk. For example if we fix k = 1 as we increase the value ofL = √ N we are moving toward lower values ofθ and the coefficient begins to rise. This phenomenon can be seen in Fig. 2 and continues for the data points not shown in the plot at smallerk/L. A similar but hierarchically less pronounced increase is observed for data points approachingk/L = 1 2 , 1 3 , 1 4 , 1 5 , 2 5 . The increase for small values ofk/L flattens out if we multiply the coefficient by (k/L) γ with γ in the range 2.5 to 3 (see Fig. 4). Since the value of F 2T has been multiplied byL 4 , the conclusion is that even if we fixk = 1 the function F 2T goes to zero in the large N limit although at a slower rate 1/L 4−γ . This question recalls the problems observed in the non-perturbative simulations of the TEK model atk = 1. Simulations at intermediate values of the coupling show a breakdown of center symmetry, which disappears when taking the large N limit at fixedθ [47]. At fixed order in perturbation theory the breaking does not take place, but the size of the corrections also points towards the benefits of keepingθ within a reasonable range. A similar analysis can be carried with respect to the potential divergences of V eff F 2T approaching the main harmonics p/q of an analogous musical scale (small q). Again the rise flattens when multiplying by (k/L − p/q) γ with the same γ as before. Once more, one concludes from this that F 2T does not diverge when taking a sequencek/L of values converging to p/q. Curiously, if one takesk = p andL = q with small q, the results have a similar size to the rest. For example for N = 4 andk = 1 or N = 9 andk = 1, which correspond tok/L = 1/2, 1/3, the values one gets for various L are not particularly large. Now we proceed to analyze the situation for larger R × R loops. The results are consistent with F 2T decreasing with the effective volume, but thek/L scaling of V eff F 2T is much less clear. One possible explanation is that as we increase R the asymptotic regime is achieved at larger values of L. As an example we we display in Fig. 5 the case of SU(4). It is clear that for R = 4 only for the largest sizes one can observe a linear approach to zero. Another aspect is also clearly illustrated by this figure: the growth of F 2T with R. Re-scaling the data by 1/R 4 we can put all data in the same plot. To see if this phenomenon extends to all values of L,L andk we studied V eff F 2T /R 4 . At fixed value of L we averaged this quantity over all values ofL andk such thatk/L ∈ [0.15, 0.45]. The filter eliminates the growth effects reported earlier for the plaquette. The final average is presented in Fig. 6 as a function of L. The results for different R are slightly displaced for visualization purposes. The error bar is the dispersion of the set of averaged values. The main conclusion is that all the values of R and L give results which are roughly of the same size of order 0.01. This is non-trivial given that the average values of N 2 F 2T have been multiplied by L 4 /R 4 ranging from 1296 to 0.0039. This leaves no doubt that the function F 2T for larger loops also goes to zero when either L or N go to infinity. Concerning the behaviour of the loops in the planes belonging to the S 2 set (02 and 13), the results are qualitatively the same as those reported previously, but the corresponding F 2T function is typically a factor two or three smaller than the one for the planes in the S 1 set.
Finally we should comment about the imaginary part of the Wilson loop coefficient described by the function G 2T . The conclusion is that this function also vanishes for either large volumes or large N . In the first case the values drop at a faster rate compatible with L −6 . Another difference, is that while the real part F 2T is typically positive, the imaginary part alternates in sign for the different values of L,L and k. The sign flips seem to coincide with the points wherek/L approaches a rational fraction p/q with small denominator, where as we saw earlier F 2T had peaks. These points coincide with those corresponding to small values of k. A good way to display our results is that given in Fig. 7. We multiplied the G 2T function by the combination of factors given below: The result for all values of L,L andk lies within a band stretching from -0.01 to 0.01. The aforementioned L dependence can be deduced from this plot. Notice that if we take the large N limit keeping k/L fixed, G 2T goes to zero as 1/N 2 . However, of one keeps k fixed and takesL = √ N to infinity, the function is going to zero at a slower rate 1 N k 2 . This matches our conclusion [47] driven from non-perturbative considerations that it is better to take the large N limit keeping k/L andk/L fixed and sizable. Our final discussion affects the sum of all contributions and the comparison of the numerical value of the second order coefficientsW (L, N, k) with that for infinite volume and number of colours. For the PBC case the leading corrections are associated to finite N or finite volume V , which go as 1/N 2 and 1/V (modulo logarithmic corrections). The coefficients are −(F 2 (∞)+F 1 (∞)/4) and −F 1 (∞)/8−R 4 ( 1 64 +(γ 2 +γ 2 log(L))) respectively. For the plaquette the first coefficient is −0.028 and the second one is of similar size for typical values of L. As R grows, the relative importance of the finite volume correction grows since it contains a term that goes as R 4 , instead of R 1.35 which is roughly the dependence of the 1/N 2 coefficient.
In the twisted case the finite volume and N corrections are blended. Keeping only the leading terms in 1/N 2 one gets Notice that for large L the first term goes to zero while the second term converges to the finite 1/N 2 correction of the periodic case. In the opposite extreme for the TEK model (L = 1) the first and second terms combine to give the finite volume correction of the periodic case. Let us now consider the general case. We may ask ourselves what configuration gives the smallest corrections at fixed number of degrees of freedom V eff ≡ V N 2 . It is clear that the first term does not depend on how we split these degrees of freedom onto spatial and colour ones. According to the analysis presented in the previous paragraphs, F 2T has a similar structure to the first term ∝ ξ/(N 2 L 4 ) with a coefficient ξ which varies slightly with the L-L splitting and the value ofk. On the contrary the second term gets smaller for larger N . In conclusion, the smallest corrections are obtained with the fully reduced TEK model, although the benefits decrease as R grows. To give a quantitative idea of the implications we see that for the plaquette expectation value the correction is 4  6 Additional considerations

Comparison with numerical simulations
Apart from the perturbative calculation we also measured the expectation value of square Wilson loops using Monte Carlo simulations. The purpose is to determine the region of values of λ = 1/b, for which this truncated perturbative expansion is a good approximation. Our methodology is based upon the auxiliary field method [67] followed by overrelaxation [68]. The numerical values of the perturbative coefficients for the twisted case are very close to those of infinite N and volume. To notice a significant effect one has to consider small values of N , L and large values of the loop size R. We first studied N = 49 with  k = 2 and measured the spatial average of the Wilson loops. To display our results instead of plotting the expectation value of the Wilson loop directly, we substract its perturbative contribution for infinite N and volume as follows: Thus, this quantity measures both the difference between the coefficients at finite and infinite N , L, as well as the effects of higher terms in the perturbative expansion. In Fig. 8 we display the result for N = 49 together with the analytic corrections to order λ, λ 2 and λ 3 . The first two come from our calculation in this paper. The latter is the result of a fit leaving the coefficientŴ 3 free. The result for the TEK model L = 1 (Fig. 8a) for the 4 × 4 loop shows that the data follow our perturbative calculation up to λ ∼ 0.15. For higher values of λ a non-zero value ofŴ (4×4) 3 ∼ 0.0195(4) is needed to match the measured value. On the other hand for L = 4 (Fig. 8b) one sees that the numerical results are unable to distinguish the first two coefficients from those of infinite N and L. The numerical value ofŴ = 0.00347 (11). The errors do not include systematics from neglecting higher orders. Unfortunately, these coefficients are not known at infinite value of N and L except for the plaquette [12,69] givingŴ (1×1) 3 = 0.000794223. We attempted a more detailed analysis in order to verify the breakdown of CP and cubic invariance induced by the twist. These effects can be seen in our calculated coefficients  Table 3: Second order coefficients of the R × R Wilson loops with R = 1 · · · 4, for the TEK model and: N = 16 (k = 1) and N = 49 (k = 2). The coefficients have been computed with quadruple precision. displayed in Table 3 for the aforementioned N = 49 k = 2 case and for N = 16 k = 1. Even for this low values of N the effects are so tiny that one needs a very high statistics study to be able to observe this breaking explicitly. For that purpose, we generated 500000 configurations of the TEK model in each case for 5 values of b (2,4,6,8 and 10). The effect is of course more pronounced the smaller the value of N and the bigger the value of R = T . We fitted the results of our Monte Carlo to a polynomial of third degree in λ = 1/b, but fixing the first two coefficients to the analytic result. This was done for the real and imaginary parts of the Wilson loops in each plane separately. The two free parameters of the fit measure the quadratic and cubic coefficients of the polynomial in λ. For the N = 49 case, the results for the quadratic piece coefficient agrees with the results of table 3. Unfortunately, the errors are of the same size as the breaking of the cubic symmetry so that this aspect could not be tested with the only exception of the imaginary part of the 4 × 4 loop. The value of this coefficient obtained for the S1 planes was 0.000385 (16), and for the S2 planes 0.000488 (36). This shows clearly both the CP and cubic invariance violation with statistical significance in agreement with table 3. In the case of the real part, although unable to show a clean plane dependence, the results were perfectly in agreement with the same table. The fitted coefficients for the S1 planes were 0.005092 (12),-0.01664(6), -0.08829 (8) and -0.20619 (12) for R = 1,2,3 and 4 respectively.
In order to see the violation of cubic invariance more neatly we also studied the N = 16 k = 1 case. Here the imaginary part (which vanishes for R = 4) shows clearly the breaking for R = 1, 2 and 3. For example for the 3 × 3 loop , the fitted coefficient is 0.00137(2) for the S1 planes and 0.00099(3) for the S2 planes. In the case of the real part there is a signal of breaking for the 4 × 4 loop, giving 0.15830(3) and 0.15778(6) for S1 and S2 respectively.

Addition of fermions in the adjoint
A very simple extension of our work is that of including fermions. There is a difficulty in including fermions in the fundamental representation since the twisted boundary conditions are singular for them. There are two ways to circunvent this problem. One is to include flavour to compensate for the boundary conditions. The other one is to allow the fermions to live in a larger lattice where they are insensitive to the boundary conditions. On the other hand there is no problem in adding fermions in the adjoint representation. There are many reasons for considering this theory interesting. One is certainly supersymmetry, but another one is the proposal done by several authors of restoring volume independence for the periodic boundary conditions case [70].
Another incentive for considering fermions is the simplicity of adding them. At the order that we are working the contribution to Wilson loop expectation values comes through a fermion loop term in the vacuum polarization, which is rather simple to add. However, the addition also induces a proliferation of options: fermion masses, number of flavours, type of lattice Dirac operator, etc. The comparison and analysis is very interesting, no doubt, but it opens up a non-trivial addition to this, already long and complex work. Hence, we opted for a mild inclusion in which we simply stick to Wilson fermions with a fixed value of the hopping parameter. The contribution of fermions to the Wilson loop amounts to the addition of a new term to the second order coefficientW TBC 2 (N, L, n µν ), which we label N f H 2 (κ, N, L, n µν ), with N f the number of adjoint flavours. Given that there is no contribution to first order there is an apparent conflict with the claim that this addition restores volume independence. However, we recall that the calculation in the case of periodic boundary conditions is not complete. We have expanded around the nontrivial holonomy ground-state and ignored the contribution of zero-modes. The addition of fermions is expected to affect the degeneracy of classical vacua which is responsible of the zero-modes.
We will now present our result for H 2 (κ, N, L, n µν ) and discuss its structure. For simplicity we will focus on the case of a symmetric box and symmetric twist in 4 dimensions. The Fourier decomposition of the adjoint fermion fields is similar to that of the gauge fields and the Feynman rules, presented in App. A, are easily derived. They lead to two extra terms in the vacuum polarization, given in App. B.1. Our expressions can be mapped to the standard ones for fundamental Wilson fermions in infinite volume [71] by performing the substitution given in Eq. (4.2) and taking into account the change in the trace normalization of the fermion representation. One of the self-energy terms is a lattice tadpole, given by Eq. (B.8). The other, Π f 2 (q) in Eq. (B.9), is the lattice analog of the fermionic contribution to the gluon self-energy. These two terms contribute at second order in λ to U (2) 2 through Eq. (3.22). They are proportional to F 2 and hence of purely non-abelian nature. Following the same strategy as in Eqs. (4.16) -(4.17), they can be decomposed in two functions in terms of which the fermionic contribution toW TBC As for the pure gluonic case, F f 2T should go to zero both in the large N and in the infinite volume limit.
We will briefly discuss below the results of the numerical evaluation of F f 2 and F f 2T for r = 1 massless Wilson fermions. Note that we can directly work with massless adjoint fermions due to the absence of zero-modes in the twisted box. Let us start by analyzing the behaviour of F f 2 at large volume. As mentioned above, this function comes from the contribution of two fermion self-energy terms. Both of them have a leading 1/L 2 correction that arises from a constant, volume-independent, term in the vacuum polarization. The structure of the correction is hence identical, modulo an overall coefficient, to the one coming from the measure. The leading 1/L 2 correction to F f 2 takes thus the form: with the same constantγ appearing in Eq. (4.20). An easy way to determine the coefficients C f i is to compute the vacuum polarization at vanishing external momentum. This is a single momentum sum whose volume expansion can be obtained following the strategy described in Appendix C. The constant, volume-independent, term is given by the infinite volume expression. In the particular case of massless r = 0 Wilson fermions, it is easy to see that it vanishes [71]. The same holds for other values of r, implying that C f 1 + C f 2 = 0.
Although not required for computing the expectation value of the Wilson loop, one can easily determine the tadpole coefficient analytically in the massless case from the infinite volume formula: with M (α) = rd − r µ cos(α µ ). The integral can be estimated numerically. In four dimensions for r = 1 we obtain C f 1 (r = 1) = −0.0612733799 (1).
With the cancellation of the leading 1/L 2 correction, the large L expansion of F f 2 in four dimensions is given by:   Notice that the coefficients of the leading fermionic and gluonic logarithmic corrections are in both cases almost independent of the loop size and opposite in sign, with the fermions counteracting as expected the gluonic contribution.
The remaining function F f 2T is very similar in structure to its gluonic counterpart but with opposite sign. It tends to zero in the same way when either N or L go to infinity. As an illustration, we plot in Fig. 9 the quantity −V eff F f 2T /4 as a function ofk/L, for the plaquette in a S 1 plane. The plot corresponds to massless r = 1 Wilson fermions. The factor 1/4 has been chosen to obtain a result comparable to the gluonic contribution. This is illustrated by displaying in the plot the pure gauge results for the L = 1 TEK model from Fig. 2. At a given value ofk/L, the two functions have the same magnitude. As a last remark, we point out that the function F f 2T for other square loops scales as R 4 like in the pure gauge case.
Although it would be interesting to explore the dependence on the fermion mass and extend this analysis to other kind of lattice fermions, this is a lengthy project that is beyond the scope of this paper and will be addressed elsewhere.

U(N) versus SU(N)
It is interesting to compare the perturbative expansion of these two groups. In the large N limit the two groups differ only by 1/N 2 corrections. In principle the U(N) group is neater as exemplified by the 't Hooft double-line notation. Our calculation was done for the SU(N) group, so that it would be interesting to know which of the 1/N 2 corrections are attributable to the restriction to this group. At leading order the result is rather simple: all 1/N 2 dependence disappears when studying U(N) instead of SU(N). Thus, the leading order coefficient is F 1 (L) for periodic boundary conditions and F 1 (L eff ) for twisted ones. - L=1, V eff F 2T gauge Figure 9: The function −V eff F f 2T (L)/4 for the plaquette in a S 1 plane is plotted as a function ofk/L for massless r = 1 Wilson fermions. For comparison we also display the pure gauge results for the TEK, L = 1 model. This is consistent with 't Hooft topological expansion which holds for U(N). All corrections in powers of 1/N 2 are associated to non-planar diagrams, which are absent at leading order.
If we proceed to next-to-leading order and focusing on the coefficient of the logarithm of the Wilson loop, one sees that the additional U(1) gluon present in U(N) only contributes to tadpole-like terms. Revising our calculation we can easily identify in what places we omitted the possible contribution of that gluon. Indeed, this was implemented by the restriction in the sum over momentum denoted by a prime. This only appeared in the terms that we called abelian: the measure contribution and the tadpole. Only the latter is affected by the addition of the U(1) mode. In summary, the second order coefficients of the logarithm of the Wilson loop for the U(N) theory are given bỹ and This matches nicely with the identification of non-planar diagrams in the U(N) theory with the only exception of the measure insertion.

Conclusions
In this paper we have studied the perturbative expansion of Wilson loops up to order λ 2 = g 4 N 2 for lattice Yang-Mills fields (with Wilson action) in a finite box with irreducible twisted boundary conditions. 5 Contrary to the case of periodic boundary conditions, this perturbative expansion at finite volume is perfectly well-defined. This is due to the absence of zero-modes. Our general presentation is valid for any irreducible twist and any dimension. The final formulas are given in terms of finite momentum sums. The effect of the different twists sits in the range over which these momentum sums run and the particular form of the momentum-dependent structure constants.
We have then studied with special focus the case of symmetric twist in a symmetric box. In particular we have analyzed the difference between the results with twist and those obtained in a simplified version of periodic boundary conditions in which the effect of zeromodes is neglected. These results depend on some common functions, F 1 (L) and F 2 (L), of the lattice size. Their large L dependence and infinite volume value have been determined for the four-dimensional case with an increased precision with respect to previous determinations. For the twisted boundary conditions case, the coefficient of the perturbative expansion to order λ 2 for each type of loop also depends on a complex function F 2T (L, N, k) + iG 2T (L, N, k). This function contains the contribution of non-planar diagrams and vanishes when either L or N go to infinity. The latter fact being a manifestation of the phenomenon of volume independence, while the former signals the independence of boundary conditions for large volumes V . These functions contain all the dependence of the result on the common flux k of the symmetric twist and, hence, are the only terms where CP and cubic symmetry violations show up. For the four-dimensional case we have evaluated the functions for a large number of values of the arguments (L, N and k) in order to determine their value and the rate of decrease to zero with either N or L. The best way to describe our findings is by plotting the values as a function ofk/L, wherek is the congruent inverse of k moduloL ≡ √ N . It turns out that for generic values of this ratio, the function decreases as one over the effective volume V eff = N 2 V . The coefficient multiplying 1/V eff grows whenk/L tends to rationals with small denominators, effectively reducing the power of N at which the functions vanish. These results are consistent with the requirement, established in Ref. [47] on the basis of non-perturbative arguments, that both this ratio as well as k/L should be kept large enough when taking the large N limit. It is interesting to realize that, although centre symmetry cannot be broken at finite N and L, our analytic calculations reinforce the interest of approaching the large N limit following our criteria.
All our results apply as well for the one-site Twisted Eguchi-Kawai reduced model (L = 1). Indeed, the best determination of the infinite N and infinite volume plaquette expectation value for a fixed finite number of degrees of freedom is obtained by using this fully reduced model. For expectation values of large loops, this advantage with respect to partial reducion (L > 1) diminishes.
We should emphasize that many of the calculations of this paper have been performed independently and using different programs by a subset of the authors. This minimizes the possibility of errors. Furthermore, we have also compared our results with Monte Carlo simulations at large values of b. A very high statistics study is necessary to verify some of the specific predictions of our calculation, such as the pattern of cubic symmetry breaking and the non-vanishing imaginary parts for each fixed plane of the loop. We found perfect agreement. Furthermore, these numerical results allow us to find out the range of values of b for which O(λ 2 ) provide a good approximation. An estimate of the O(λ 3 ) coefficients has also been obtained.
A bunch of additions have been included to make this paper as complete as possible. In particular, we have analysed the difference between the U(N) and SU(N) cases, and more importantly we have also computed the effect of including fermions in the adjoint representation. These fermions are fully compatible with twisted boundary conditions and have been subject of great interest because of their role in supersymmetry, volume independence [70,[72][73][74][75][76][77][78][79][80][81], infrared fixed points [82][83][84][85][86], etc. In relation with our perturbative results, their contribution only enters via the self-energy of the gluons and is proportional to the number of flavours N f . For the four-dimensional symmetric twist case, we studied its effect on F 2T for massless Wilson fermions. The result is similar qualitatively and quantitatively to the purely gluonic results. Given their interest a more detailed analysis using different versions of lattice fermions and different masses is well justified but falls away from the main scope of this paper. this condition the δ functions imposing momentum conservation at the vertices should be understood in strict sense and not modulo 2π.

B The vacuum polarization at O(λ 2 )
For completeness we give the expression for the vacuum polarization for the Wilson action up to order λ 2 as derived by Snippe in Ref. [38], generalized to a twisted box with an arbitrary irreducible twist.

B.1 The contribution of adjoint Wilson fermions
The fermionic contribution to the vacuum polarization is proportional to the number of fermion flavours N f . For N f = 1 adjoint Wilson fermions and up to order λ 2 , it includes two terms given by: where N σδ (q, p) = δ σδ cos 2 (2p + q) σ 2 M (p)M (p + q) + ρ sin p ρ sin(p + q) ρ (B.10) Under the substitution Eq. (4.2) and the change in trace normalization of the fermion representation, our formulas reproduce the infinite volume results by Kawai, Nakayama and Seo in Ref. [71].

C Evaluating finite volume corrections
Here we consider the evaluation of quantities of the type where the momenta p are d-dimensional vectors where each component has the form p µ = 2πmµ Lµ with integer m µ ranging from 0 to L µ − 1. The sum extends over all values of m µ except when all of them vanish simultaneously. Finally V stands for the volume, which is equal to the product of all L µ . The goal would be to obtain the large volume behaviour of I(L), and in particular the corrections to the infinite volume limit.
Traditionally in dealing with sums of our type one makes use of the Euler-MacLaurin formula. Here however our integrands are periodic (with period 2π) in each variable. To treat this type of integrals we use the following expression where Dα is the product of dα µ /2π over all directions. The next step is to sum over p, giving Plugging the expression onto the formula we get The last term is equal to F(0)/V , unless F(0) diverges. Notice that the term in the first sum corresponding to n µ = 0 for all µ gives the infinite volume limit of the expression, which we assume to be convergent. This is the leading contribution in the Euler-MacLaurin formula. Thus, we get an exact expression for the finite volume corrections to our integral as follows where the prime means the sum over all d-dimensional vectors n ∈ Z d , excluding n = 0. The argument of the exponential is a simplified form meaning µ n µ L µ α µ .

C.1 The finite volume propagator
Now we can apply our formalism to the study of the expectation values of Wilson loops. Our main integral under consideration comes from using as integrand the following expression whereD(α) = 4 µ sin 2 (α µ /2)+m 2 = 2d−2 µ cos(α µ )+m 2 . For β = 1 the corresponding sum I is nothing but the propagator of a scalar particle of mass m on a finite lattice. The mass is necessary to have a well defined value of F(0). In all the main expressions that we will use later it would be possible to take the limit m −→ 0. Working with β different from 1 allows to evaluate other expressions (like the measure contribution to the vacuum polarization given by β = 2) and can also act as a regulator. Now we can use Schwinger trick to recast the integrand as an integral as follows The factorized exponentials can then be treated as the integrand and we can apply our formalism to it. The infinite volume limit result is then given by where I l is the modified Bessel function. For β = 1 this is just the lattice propagator at infinite volume. The finite volume propagator has then the form of the sum of the propagators to all replica points l µ + n µ L µ . The lattice propagator at large distances is well-approximated by the continuum one. This follows from the asymptotic expansion of the Bessel functions at large values of x. The leading term is as follows: Thus, the leading finite volume correction to P (l, β) is given by where z µ = l µ /L µ and the function G(x, z, L) is given by where ϑ(z; iτ ) is the Jacobi theta function, whose duality relation has been used for the last equality. The first term on the right-hand side of Eq. (C.8) is just the F(0)/V subtraction. When the sizes go to infinity uniformly as follows L µ = λ µL , it is clear that the integrand of Eq. (C.8) is strongly suppressed whenever x/L 2 1. Thus, we can change variables to y = x/L 2 and restrict the integral to go from to infinity. This gives (C.10) The convergence of the integral at small y is guaranteed by the behaviour of the quantity inside parenthesis irrespective of the other factors. At large y the mass term guarantees convergence. Taking the mass to zero produces a divergence of the integral, coming from the large x behaviour of ϑ(z; i 4πx L 2 ) −→ 1. Its value is subtracted by the leading first term giving a convergent result. Hence one can combine these terms and give a result after take the limit of vanishing mass and in this combination. We have still preserved the mass dependence in the last term. It is necessary to ensure its convergence at large y whenever d ≤ 2β. This is the same range for which the massless infinite volume limit does not exist. Outside this range one can savely set m = 0 in Eq. (C.11).

C.2 Leading order of Wilson loop
Now we can apply the formalism to the expectation value of Wilson loops at lowest order. The integrand is given by F(α) = sin 2 (α 0 T /2) sin 2 (α 1 R/2) sin 2 (α 0 /2) (D(α)) β + ( R↔T 1 ↔ 0) with β = 1. It is possible to write the expression as follows which in the presence of a mass term vanishes at α = 0. Notice that we take the loop in the 0 − 1 plane with size T and R respectively.
We can relate the calculation to the previous one of the P (l; β, L) We introduce the displacement operator δ 1 which adds 1 to l 1 . We call δ −1 1 the inverse operator which displaces by −1. In an analogous fashion δ 0 displaces l 0 . With this notation the finite volume correction to the leading contribution to the wilson loop is given by In order to obtain the leading result it is interesting to consider the limit in which z 1 = R/L 1 is treated as a small quantity. The operator can then be expanded as where ∂ 1 is the derivative with respect to z 1 (treated as a continuum variable). If we now apply the same procedure to the operator along the time direction, we get We are now in position to compute the leading correction to the Wilson loop. All we have to is to apply the operator to our previous expression P (l, β; L) and then set z = 0. We can make use of the result   +0 ↔ 1 (C.12) in which the mass has been set to zero, assuming that β − 1 < d/2. In the particular case in which all lengths are equal L µ =L, there is considerable simplification since the expression inside the parenthesis becomes a total derivative For the Wilson loop (β = 1) the result is just given by the value of the function at the limits ( -1 at y = 0) giving δW (1) = − T 2 R 2 2dV (C.14) For the measure term (β = 2) the formula for the symmetric case can be obtained by integration by parts and gives F mes (L) − F mes (∞) = − δW (2) 12 = − T 2 R 2 24dL d−2 ∞ 0 dy ϑ d (0; i4πy) − 1 − 1 (4πy) d/2 (C.15) The two dimensional case is the only physical case for which the previous expression is not valid. This is so because there is no massless infinite volume limit. To deal with this case one has to go back to the expression including a non-zero mass. The divergence comes from the last term which when integrated from y = 1 gives the incomplete gamma function Γ(0, m 2 L 2 ) = − log(m 2 ) − log(L 2 ) − γ + . . .. The mass singularity cancels with that coming from the infinite volume limit quantity and shows that F mes (L) = − R 2 T 2 96π log(L) + constant + . . . (C. 16) which is used in the text.

D Non-abelian contributions
In this section we develop the methodology to study L and N dependence of expressions of the form where the momenta p = p s + p c , with p s and p c are the spatial and color momenta respectively. The integrand is an unspecified function A(p, q), which is assumed to be periodic of period 2π in each of the arguments and regular everywhere. Finally F (p c , q c , −p c − q c ) is the characteristic structure constant in the colour momentum basis. We recall that with the hermiticity convention that we are using one has F 2 (p c , q c , −p c − q c ) = 1 N (1 − cos(2θ(p c , q c ))) = 1 N (1 − cos(Φ(p c , q c ) − Φ(q c , p c ))) (D. 2) It is convenient to consider the two terms separately as we did in the analysis of the results performed in section 4. The first part contributed to all functions labelled N A. They are called like that because they arise from terms involving the F 2 , which are structure constants of the non-abelian group. In this case there is no colour factor and the corresponding sum will be given by The part containing the cosine was labelled NP, standing for non-planar, since this twist dependent factors only occurs in the non-planar part of the diagram. The resulting expression is I N P (L) = − 1 V 2 eff pc,qc∈Λmom/Λ L ps∈Λ L qs∈Λ L A(p, q) cos(2θ(p c , q c )) (D.4) To start with, let us consider the infinite volume limit of the non-planar part. We can use the Euler MacLaurin result, or the analysis performed in the previous appendix to conclude that it is given by which is just the well-known value of the quadratic Casimir in the adjoint representation. Using Eq. (D.9) one can obtain the infinite volume limit of the non-planar part, given by which is just the infinite volume limit of the original expression but omitting the colour degrees of freedom.
One can go beyond this result and try to evaluate the finite volume corrections to this non-planar part. We can use the formalism introduced in the previous appendix to replace all sums over space momenta by integrals. Finally, we get The previously obtained infinite volume result corresponds to takingm µ =ñ µ = 0 for all µ. Excluding this value from the sum we get the finite volume correction. Let us process our result a bit more by realizing that the cosine only depends on the arguments modulo N. Hence, one can split the integers as followsm µ = l µ + Nn µ . Indeed for the symmetric twist cases the argument applies withL replacing N . Thus, we can rewrite The formula is valid for any twist if we identifyL with N and n µν with the twist tensor. However, for the four dimensional symmetric twist it is more convenient to takeL = √ N and n µν = n µν /L.
The function H is an oscillatory function with periods proportional to 1/L. Let us now restrict ourselves to the symmetric twist case in a symmetric box of size L in both 2 and 4 dimensions. In that case the tensor n µν = k µν where k is an integer coprime witĥ L and µν is an invertible antisymmetric matrix. We might redefinel µ = µν l ν . Due to the invertibility, the range over whichl µ runs coincides with that of l µ . We can also change variables from β toβ given byβ