The twisted gradient flow coupling at one loop

We compute the one-loop running of the SU(N) ’t Hooft coupof twelve-flavor SU(3) gauge theoryling in a finite volume gradient flow scheme using twisted boundary conditions. The coupling is defined in terms of the energy density of the gradient flow fields at a scale l˜\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{l} $$\end{document} given by an adequate combination of the torus size and the rank of the gauge group, and is computed in the continuum using dimensional regularization. We present the strategy to regulate the divergences for a generic twist tensor, and determine the matching to the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} scheme at one-loop order. For the particular case in which the twist tensor is non-trivial in a single plane, we evaluate the matching coefficient numerically and determine the ratio of Λ parameters between the two schemes. We analyze the N dependence of the results and the possible implications for non-commutative gauge theories and volume independence.


JHEP03(2019)200 1 Introduction
In recent years, the continuous smoothing procedure known as the gradient flow [1][2][3] has received considerable attention. One of its most common applications has been, in combination with finite-size scaling techniques, the determination of the non-perturbative scale dependence of the gauge coupling constant. Examples of the usefulness of this approach range from precise, non-perturbative determinations of the QCD coupling constant and Λ parameter [4][5][6][7][8] to the study of Yang-Mills theories with near conformal behavior [9][10][11], or with large number of colors [12]. Several coupling renormalization schemes based on gradient flow techniques have been proposed to that end [3,[13][14][15][16]. These schemes can be related through a perturbative calculation to more traditional ones such as the MS scheme, a step often required to make contact with experiment. However, and despite their importance, perturbative calculations in such a set-up are scarce. In infinite volume, the matching has been determined up to both next-to-leading order (NLO) [3] and next-to-next-to-leading order (NNLO) [17], while for finite volume it has been done up to NNLO in the Schrödinger functional scheme using numerical stochastic perturbation theory [18,19].
In the scope of this work, we will focus on a particular gradient flow finite volume scheme for SU(N ) pure gauge theory, introduced by A.Ramos in ref. [16]. We will be presenting results for the matching at NLO of this scheme to the MS one, determining the coupling scale in terms of the size of a 4-dimensional torus endowed with twisted boundary conditions (TBC) [20]. From the point of view of perturbative calculations, TBC have an enormous advantage over periodic ones (PBC), as using TBC turns the set of zero-action solutions into a discrete one, and avoids the quartic nature of the fluctuations around A µ = 0 present with PBC [21]. The usefulness of TBC for perturbation theory was first formulated in the context of volume reduction in large N Yang-Mills theory [22,23], and was then extended to various other contexts at finite and large N [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. Despite these advantages, as we will show along this work, perturbative calculations in the twisted gradient flow scheme remain challenging, although so far an analogous perturbative calculation in the case of periodic boundary conditions has not been obtained.
Our interest in this calculation goes beyond the particular applicability of the results, and connects with theoretical ideas related to the concept of volume independence in gauge theories. An essential ingredient of the construction has to do with the dependence of the coupling on the number of colors. We will follow the finite-volume prescription adopted in [12], and set the SU(N ) running coupling scale to be proportional to an effective sizel combining the torus size and the number of colors. In perturbation theory, this effective scale is expected to jointly capture the dependence of the coupling on the volume and on N , once an angular variable depending on the choice of twist is fixed. For the purposes of this paper, we will keep N finite but use our results to analyze the dependence of the coupling on N , at large values of N . Setting the energy scale of the renormalized 't Hooft coupling to be µ = 1/(cl), we will consider two different types of N → ∞ limits, a thermodynamic limit in which the effective size is sent to infinity as c goes to zero while µ is kept constant, and a second one, called singular in [39], in which N is sent to infinity as the torus size is

JHEP03(2019)200
shrunk to zero at constant c in such a way as to keepl fixed [12,34,37,40]. 1 In ref. [12] this last limit was used to compute the SU(∞) running coupling through a step scaling procedure in which the step size was modified via changes in the rank of the gauge group.
This particular prescription for scale setting was inspired by the idea of volume reduction in lattice gauge theories. Originally formulated by Eguchi and Kawai [45], volume reduction states that in the (thermodynamic) large N limit, SU(N ) theory becomes independent of the physical size of the torus. Proof for this statement relies on the independence of the large N Schwinger-Dyson equations from lattice volume, which in turn requires center symmetry to be preserved. As the symmetry was shown not to hold with PBC [46], several alternative proposals were formulated [42,[46][47][48][49], one of which was the use of twisted boundary conditions [22,23], which has proven very successful provided the twist tensor is judiciously chosen [35,36,[50][51][52][53][54][55]. The idea of volume reduction with TBC was extended to the continuum theory in ref. [56], constituting the first formulation of the Feynman rules of a non-commutative Yang-Mills theory [57]. The singular limit in ref. [39] was first formulated precisely in the context of such non-commutative theories, with the effective torus sizel arising in a natural way through the Morita duality. Thus, ordinary gauge theories on a twisted torus are related to non-commutative theories with a rational value of the dimensionless non-commutativity parameter, the effective size corresponding to the one of the non-commutative torus. In 2+1 dimensions, several recent works [37,40] have analyzed the possibility of defining non-commutative gauge theories at irrational values of the non-commutativity parameter as the limit of a sequence of ordinary twisted gauge theories of an increasing number of colors. These works have shown that, if one wishes to avoid tachyonic instabilities [58], such a construction can only be achieved for an uncountable, zero-measure set of values of the non-commutativity parameter. We will, in light of these results, analyze the behavior of the coupling in the singular large N limit.
In this paper, we present a perturbative calculation in the continuum of the running 't Hooft coupling constant at NLO. The layout of the paper is as follows: in section 2, we introduced the twisted gradient flow (TGF) scheme, presenting the gradient flow observable used to define the running coupling (i.e. the energy density evaluated at a positive flow time proportional to an effective sizel), along with some specifics about the implementation of TBC in our setting. In particular, we detailed the orthogonal twist used throughout the paper. Section 3 then presents the perturbative expansion, along with the regularization and renormalization schemes. This is the longest, most technical section, and contains a fair bit of algebraic manipulation. The calculation is analogous to the one performed by Lüscher in infinite volume [3], though many particularities to the twisted finite volume scheme appear. We will simply mention that it contains the expansion of the observable in powers of the coupling, a reformulation of the NLO contribution as the sum of several integrals, the identification of the divergent terms entering the calculation, and a procedure to regularize them by relating them to infinite volume expressions that can be evaluated in dimensional regularization. Expressions for the observable at LO and NLO are provided JHEP03(2019)200 within the section, but any reader interested in the final expression for the matching at one-loop order of the TGF coupling to the MS scheme may skip directly to section 4, which contains both the matching to the MS scheme and the ratio of Λ parameters (which need to be computed numerically). Results for the case of a 2-dimensional non trivial twist and several SU(N ) groups are presented in subsection 4.2. Section 5 discusses the dependence of the coupling on the number of colors, following similar arguments to those in [36,37]. A summary of results is presented in section 6. Many technicalities were moved for clarity to appendices A-D, including details on the algorithms used to compute the Λ parameter.

The twisted gradient flow coupling
One of the applications of the gradient flow method has been the computation of the Yang-Mills running coupling, using the energy density E(t) of the gradient flow field as the defining observable. At positive, non-zero flow time t, t 2 E(t) is a renormalized quantity and, at leading order in perturbation theory, is proportional to the MS coupling at a scale µ = 1/ √ 8t, which leads to a natural definition of a renormalized coupling constant [3]. In this work, we focus on a particular gradient flow scheme that makes use of finite size scaling on a torus with TBC. As discussed in the introduction, our set-up is based on the one introduced by A. Ramos in ref. [16], but differs slightly from it for reasons that will become clear in what follows.

The definition of the coupling
The gradient flow is based on the introduction of a parameter t, known as flow time, in such a way as to define a t-dependent gauge field B µ (x, t) matching the Yang-Mills one A µ (x) at t = 0. As flow time passes, this gauge field is smeared down towards the minimum action solutions, its evolution driven by the so-called flow equations: where D µ and G µν respectively stand for the covariant derivative and field strength tensor of the flow fields: This scheme is particularly useful, as observables built from the expectation values of products of B fields at positive flow time have been shown to be renormalized quantities [59]. The renormalized gradient flow coupling can then be defined in terms of the energy density of the flowed field: In infinite volume, this quantity can be used to define a renormalized SU(N ) 't Hooft coupling at an energy scale µ, given by [3]:

JHEP03(2019)200
where d stands for the number of space-time dimensions of the theory. Finite volume gradient flow schemes [12][13][14][15][16] use a formulation in which the gauge theory is defined on a finite torus instead, with each scheme differing in specific details such as, for instance, different boundary conditions. The most common choice is to use a symmetric torus, with all directions of equal length l, while setting the scale for the renormalized coupling in terms of l by fixing µ = 1/(cl), with c an arbitrary constant. Each specific choice of c, always taken to be smaller than unity, is an intrinsic part of the definition of the scheme. The SU(N ) TGF coupling used in this paper is inspired by the finite-volume schemes proposed in refs. [12,16,60]. We will leave the specifics of the scheme for the next subsection, but mention that our SU(N ) gauge theory will be defined on an 4-torus with TBC [20], and such that the torus has a period l in d t (twisted) directions, andl = lN 2/dt in the remaining 4 −d t ones, with d t being either two or four. The reasons behind our choice of an asymmetric torus will become clear in what follows. In this scheme, the twisted gradient flow 't Hooft coupling is defined by [12]: where F (c) is a constant defined in such a way as to have λ TGF (l) = λ 0 + O(λ 2 0 ), in terms of the bare 't Hooft coupling λ 0 .

The choice of boundary conditions and torus size
In this subsection, we will discuss the particular definition of the TGF scheme used in this paper. The main idea is to have a perturbative set-up that is as symmetric as possible [12,60]. To achieve this, we will look at the quantization of momenta in our particular setting, and select the torus size accordingly. We will begin with a generic discussion of the quantization of momenta in the presence of TBC, leading to the introduction ofl as the relevant length scale.
Let us start by considering a SU(N ) gauge theory defined on a d-dimensional torus of length l µ in each direction, and focus our analysis in the specific case of four dimensions for a gauge potential that satisfies 't Hooft TBC [20]. We will work with an orthogonal twist, for which the gauge potential can be fixed to be periodic in each direction up to a constant gauge transformation: where Γ ν are four SU(N ) matrices known as twist eaters which satisfy: with Z µν an element of the center of the gauge group, written in terms of an antisymmetric tensor of integers n µν as: The twist tensor n µν is preserved under gauge transformations, and uniquely characterizes the boundary conditions. It is said to be orthogonal when κ(n) = µνρσ n µν n ρσ /8 = 0

JHEP03(2019)200
(mod N ). Among such tensors, we will focus only in the so-called irreducible twist tensors, which are the subset for which the only matrices that commute with all Γ µ are the ones proportional to the identity in SU(N ). Irreducible twist tensors have been known to be advantageous for perturbative calculations, as the class of gauge-inequivalent zero-action solutions is discrete [61,62] and zero-modes are eliminated, making computations in perturbation theory much easier. A detailed discussion of the conditions under which a twist is irreducible can be found in [63]. For the scope of this work, we will focus on two types of irreducible twist tensors (detailed below), which are non-trivial in either a single plane or in all of them. For the sake of clarity in the description, we will use gauge freedom to impose strict periodicity for the gauge potential in all directions except for a number d t of them, dubbed "twisted directions", taken to be either two or four, though the specific form of the twist matrices is irrelevant as long as eq. (2.8) is satisfied. We will write our orthogonal twist tensor in the form: where l g = N 2/dt depends on both the number of colors and the number of twisted directions, and k and l g are two coprime integers that guarantee that the irreducibility condition is satisfied. The choice to have a non-trivial twist in only the (0,1) plane is made by setting d t = 2 and 01 = − 10 = 1, and by choosing µν = 0 in any other plane, whereas to twist all planes non-trivially one must instead take d t = 4, and set µν to be antisymmetric and equal to 1 whenever µ < ν. With this choice: A non-trivial twist, such as the one above, will affect the quantization of momenta in the finite box.
The solution to the boundary conditions on such twisted tori in the continuum is well known [23], as one can see for instance in [36] for the general treatment when the torus is discretized on a lattice, or in [34] for an example in 2+1 dimensions in continuum perturbation theory. We will, in what is left of this subsection, recall some known results necessary to implement perturbation theory with TBC. We start by defining: with s µ (q) ∈ Z. Provided k and l g are coprime integers, there will be N 2 independent SU(N ) matrices of this type, of which the only non-traceless one is the one proportional to the identity matrix, i.e. the one for which s µ (q) = 0 (mod l g ) in all twisted directions. Excluding it, the remaining N 2 − 1 matrices constitute a basis for the SU(N ) Lie algebra. IfΓ(q) satisfies: 13) with no summation over ν implied, the boundary conditions in (2.7) are trivially implemented through the Fourier expansion: 14)

JHEP03(2019)200
where V ≡ µ l µ , and the prime in the sum denotes the exclusion of the momenta for whichΓ(q) ∝ I. In the periodic directions, for which Γ ν ∝ I, the momenta are as usual quantized in units of 2π/l ν . This is however not the case for the twisted directions, where a solution is provided by: wherek and˜ µν are given by: The momentum along the twisted µ directions is thus quantized in units of 2π/l µ , with l µ ≡ l µ l g . For this choice of s µ , the group structure constants in theΓ(q) basis become momentum dependent and are given by: and The tracelessness of theΓ(q) matrices thus forbids momenta such thatl µ q µ = 0 (mod 2πl g ) in all twisted directions, and so in particular it forbids zero momentum in the twisted box.
The previous analysis implies that momentum is quantized differently in periodic and twisted directions: it is quantized in terms of the inverse torus size for the former, and in terms of an effective size combining the torus period and the number of colors of the gauge group,l µ = l µ l g = l µ N 2/dt , for the latter. This observation has led us to a specific choice of torus size to define the TGF coupling in eq. (2.6), picked in such a way as to impose the same momentum quantization in all directions. When d t = 2, this will be achieved by considering an asymmetric torus of length l in the twisted directions andl in the periodic ones, whereas for d t = 4 we will instead pick a symmetric 4-torus of period l in all directions. This way, all momenta will always be quantized in units of 2π/l, and we will use this effective sizel as the renormalization scale for the running coupling.

Perturbative expansion
The procedure to determine the perturbative expansion of the coupling follows closely the one developed by Lüscher in infinite volume in ref. [3]. The main difference arises from the quantization of momentum on the torus, as momentum integrals become sums over an infinite set of discrete momenta, and from the change in the group structure constants due to the different choice of SU(N ) Lie algebra basis. Divergent momentum sums, however, can still be treated via dimensional regularization -see for instance [64] -in a way that will be detailed in this section.

Perturbative expansion of the energy density
As a first step towards obtaining the perturbative expansion of the observable, we will fix the gauge in such a way that the following periodicity conditions are satisfied: where Γ µ satisfies eq. (2.8), and with a twist tensor of the form shown in eq. (2.11). This restricts the set of allowed gauge transformations Ω(x) down to those preserving the form of the twist matrices, i.e. those satisfying: These boundary conditions are implemented through the Fourier expansion of the gauge field given in eq. (2.14). In the specific case of the asymmetric torus that we are considering, the torus volume is given by V = l dtl d−dt , and momenta in all directions are quantized in terms of the effective sizel. As we recall, the prime in the sum in eq. (2.14) denotes the exclusion of all momenta for whichlq µ = 0 (mod 2πl g ) in all twisted directions, which in particular excludes zero modes. With this, we may begin the perturbative expansion, which we perform around the A µ = 0, zero-action solution. We start in d = 4−2 dimensions by scaling the original gauge potential with the bare coupling, A µ (x) → g 0 A µ (x). The full Feynman rules in momentum space, given in the Feynman gauge and derived using the boundary condition-preserving Fourier representation mentioned in the previous section, can be found in appendix A.
It will be convenient to henceforth use a set of modified flow equations:

5)
ξ being a gauge parameter to be set to unity. At fixed t, the field derived from this modified flow equation can be related to the solution of the original one by a gauge transformation [3], and hence the modification does not affect gauge invariant observables such as the one we are considering in this paper. It can be shown that the corresponding (flow-time dependent) gauge transformation preserves the boundary conditions (3.1) at any given flow time [16]. These modified flow equations can be solved order by order in g 0 by expanding the flow field in powers of the coupling: The flow field satisfies the same boundary conditions as the original gauge potential and can be Fourier expanded, at any given order, in the same way:

JHEP03(2019)200
The expansion of the energy density in powers of g 0 can now be obtained by expanding the fields in E(t) directly. Dropping for clarity the arguments of the fields in position space, B . The corresponding expression in momentum space, however, is specific to the TGF set-up. In particular, the SU(N ) structure constants f abc appearing in infinite volume are replaced by the momentum dependent functions F (p, q, r) appearing in the commutation relations of theΓ(q) matrices -see eqs. (2.17), (2.18). For the sake of completeness we give below the seven different terms contributing to the expectation value of E(t) arising at order g 4 0 , with an additional 1/N normalization factor added for later convenience. Each term can be identified with one of the lines in eq. (3.8): 14)

JHEP03(2019)200
(3.15) The shorthand notation p i in the δ functions was used to denote the sum over all present momenta for each term. The E 0 term will turn out to be a combination of a leading O g 2 0 term and an O g 4 0 correction, whereas all other terms will turn out to be O g 4 0 . Then, the next step is to relate the flow fields to the actual gauge fields A µ (x), for which we will need to obtain an order-by-order solution to the flow equations.

Solving the flow equations in the TGF scheme
Let us consider the flow equation (3.5) with the gauge parameter ξ set to unity. This was already solved by Lüscher for the infinite volume case [3], but the results in finite volume are slightly different. Expanding the fields in perturbation theory, as in eq. (3.6), and dropping for clarity of notation the arguments of the fields in position space, the equations to solve order by order are of the form: The first three orders will be enough to obtain the observable at order O g 4 0 : We may define a momentum space version of R (i) µ : under which the R µ terms read:

JHEP03(2019)200
In terms ofB µ (q, t), the flow equation in momentum space becomes: whose solution is immediate at first order: and which can be solved for the next two orders by directly integrating R (i) µ : Higher order terms, while increasingly tedious, can be obtained through the same iterative procedure. From these expressions, and using the Feynman rules from appendix A, we derived the expressions of the contributions from eqs. (3.9)-(3.15) in terms of sums over momenta. Introducing for the sake of readability the symbol: and after quite a bit of algebra, we ended up with:

JHEP03(2019)200
where the bare coupling and the volume have been replaced by the bare 't Hooft coupling and the effective length, and we have defined an auxiliary momentum p = q + r. The primes from the sums in the O(λ 2 0 ) terms have been discarded, as the F 2 (q, r, −q − r) factors automatically vanish for such momenta.

The energy density at LO
As we recall, our aim in this paper was to obtain a perturbative expansion of the observable E(t)/N at NLO, which in powers of the bare 't Hooft coupling can be parametrized as: We will begin by deriving the leading order term from the formulas in the previous subsection, given by: It will be convenient to introduce a few auxiliary variables and functions: in terms of which we may write: This A function can be expressed in terms of Jacobi theta functions θ 3 as: where we used Poisson resummation to rewrite the theta functions: The leading order infinite volume expression is retrieved in the c → 0,l → ∞ limit, taken in such a way as to keep cl fixed. In that limit: leading to: 45) in agreement with the results in ref. [3].

The energy density at NLO
As for the subleading O(λ 2 0 ) term coming from eqs. (3.29)-(3.35), we found that, after a fair bit of algebra, it can be expressed in terms of a handful of integrals. By rewriting the momenta in denominators as exponents using Schwinger's parametrization, and the momenta in numerators as derivatives with respect to the flow time variables, we were able to recast the expression for the energy density at NLO as: where the I i are twelve relatively simple integrals to be detailed later on. As the computations and manipulations are rather long and tedious, we will illustrate the procedure using one of the simplest contributions, E 4 in eq. (3.33), and show the remaining E i contributions in terms of the basic integrals in appendix B. The E 4 contribution is given by: Using Schwinger parametrization to lift the momenta from the denominator, we then defined an integral: so as to rewrite: The presence of the structure constants in each E i will let us formulate the integrands in terms of Siegel theta functions. We may indeed rewrite N F 2 as: a substitution under which a generic integrand of the form: becomes: where we rescaled the variables s ≡ 4πl −2s , u ≡ 4πl −2ũ , v ≡ 4πl −2ṽ , and where we used the quantization of momenta in the twisted finite box to rewrite q and r in terms of integers.

JHEP03(2019)200
The connection to Siegel theta functions becomes clear by introducing the function In this expression I d denotes the d×d identity matrix and the sum over M denotes the sum over the corresponding integers m µ , n ν , regrouped into a 2d-dimensional column vector.
Recalling the definition of the Siegel theta functions: this matricial expression takes the form: (3.56) Using this notation, the integral entering E 4 reads: With this, only one last bit of manipulation is left in order to have the integrals ready for the calculation of the energy density at NLO. In terms of the variables t andĉ defined in eq. (3.38) and rescaling z appropriately we have: Introducing an auxiliary function Φ(s, u, v,θ) that incorporates the normalization factor in front of the integral: Φ(s, u, v,θ) = N G(ĉs,ĉu,ĉv,θ), (3.59) with: we rewrote the integrals in a fairly basic form allowing us to evaluate them numerically. For instance, for the integral in E 4 : A similar procedure can be followed for all the terms contributing to the energy density at NLO, leading to the result in eq. (3.46) for E (1) (t), where the twelve intervening integrals are:

Structure of UV divergences
As some of the integrals defined in the previous subsection are UV divergent in 4 dimensions, in this subsection we will discuss how to parametrize their asymptotic behavior. We will show that, in all of our cases, the divergent contributions can be expressed in terms of an infinite volume integral that can be regularized through analytic continuation in d.
The relation to the existing infinite volume calculation from ref. [3] will be presented in section 3.3. The UV singularities are tied to the structure of the Siegel theta functions entering the definition of the Φ function: The real part of the matrix A(ĉs,ĉu,ĉv,θ), obtained by settingθ = 0 in eq. (3.54), is a positive definite symmetric matrix as long as det A(ĉs,ĉu,ĉv,θ = 0) > 0, i.e. when (su−v 2 ) > 0, which ensures that the series defining the theta function converges uniformly. It will be useful to define a new quantity: which is always positive definite in our integration ranges, and hence the determinant will be positive definite as well, except at the points for which u = 0. 2

JHEP03(2019)200
The analysis of the asymptotic behavior of the integrals is much clearer once we apply Poisson resummation as in eq. (3.43) to each component of n in the definition of Θ: (3.76) Wheneverθ˜ m / ∈ Z d , the corresponding term will be asymptotically finite at u = 0. However, in the case in which we have a vector of integers, we will be able to remove theθ dependence by shifting n, thus leaving the asymptotic behavior to be driven by the shifted n = 0 terms. Such terms go, as we approach u → 0, as: This observation allows us to isolate the asymptotic divergence by identifying the cases for whichθ˜ m ∈ Z d . The first case in which this occurs is wheneverθ ≡k/l g = 0, for any value of m. For nonzeroθ, it will happen whenever˜ m = 0 (mod l g ). Since the vector m has nonvanishing components only along the twisted directions, this will be the case whenever m µ = 0 (mod l g ) simultaneously for all twisted directions. The terms responsible for the UV divergences at u = 0 have therefore been identified, and come in two categories: • Forθ = 0, terms with n = 0 and any value of m.
With this, we may begin the discussion on how the divergent integrals can be regularized.

Regularization
The regularization strategy will be based on splitting each integral into the sum of a finite piece that can be directly evaluated at d = 4 and integrated numerically, and an asymptotic term to be handled analytically using dimensional regularization. The way to implement such a strategy will be discussed in this subsection. We start by introducing a function H(s, u, v,θ) given by: with the usual meaning for the prime in the sum over m. The Φ function entering the integrals can then be rewritten as: which is quite advantageous, as the explicit exclusion of the momenta m proportional to l g from the sum automatically makes the term inθ = 0 finite at u = 0. All UV divergences at d = 4 thus come, in this parametrization, from the H(s, u, v, 0) term, and are of the form: with α defined as in eq. (3.75), and, as we recall, positive definite everywhere in the integrals. Hence, the sum over m is convergent, and the leading asymptotic behavior at u = 0 is controlled by the u −d/2 factor (times the additional powers of u appearing in the integrand prefactor). It will be useful to write the function Φ (0) in terms of the function A(x) from eq. (3.39): For reasons that will become clear later, we define: in terms of which we may rewrite: This formulation will be useful to analyze the asymptotic UV behavior of the integrals resulting from replacing the original function Φ in the integrand by Φ (0) . Before discussing the general treatment, we will deal with I 1 as a representative example. In this case the integral diverges at x = 0 as u = 2x. The piece containing the divergence thus reads: which, substituting in the expression for Φ ∞ , yields: The asymptotic behavior at small x can then be obtained by expanding A (2ĉt −ĉx/2) around x = 0. The integrand of the leading term goes as x 1−d/2 , whereas the next to leading term is convergent in d = 4. Hence, the integral will behave asymptotically as: Notice that in this expression the entire momentum dependence has been factorized into the normalization constant A(2ĉt ), which happens to be the same factor that appeared at leading order -see eq. (3.40). The integral I 1 (Φ ∞ , t ) can then be evaluated in dimensional regularization with d = 4 − 2 , leading to: The asymptotic expansion of all other integrals (except for I 9 , which we will address separately) is obtained in the same way: we expand the function A(ĉα) appearing in the definition of Φ (0) around u = 0, retain the leading term, and then use it to define:

JHEP03(2019)200
Remarkably, the integrals I i (Φ ∞ , t ) match the ones appearing in the infinite volume calculation (up to a factor depending on N ), which we will present in section 3.3.
We are now in a position to summarize, still keeping I 9 aside, the regularization strategy. The idea is to decompose the finite volume integrals into two pieces, one that is finite in four dimensions: and another one, shown in eq. (3.88) above, that requires analytic continuation to four dimensions and is proportional to each corresponding infinite volume integral. The ultraviolet divergences of the original integral are contained in this last piece, and appear as poles in 1/(d − 4), though only I 1 , I 4 , I 5 and I 7 turned out to have such 1/ poles. As for the strategy to regularize I 9 , some modifications, described in detail in appendix C, are required. The initial integral is decomposed as: with the Heaviside step function θ restricting the interval of integration over z. The first term on the right hand side is finite in four dimensions, while the other two have to be analytically continued to d = 4. Denoting I reg 9 these analytic continuations, we end up with: And therefore: and:

Infinite volume limit
The expression of the energy density in infinite volume can be easily retrieved (see [36]) by making the following substitutions in eqs. (3.28)-(3.35): The resulting expressions for the contributions to the energy density, after integrating over the d-dimensional momenta, can once again be rewritten in terms of twelve basic integrals,

JHEP03(2019)200
much like what happened in the finite volume case. We will first present the case of E 4 as an illustrative example, and then present the results for the general case. The infinite volume expression for E 4 is obtained making the substitutions from eqs. (3.95) and (3.96) in eq. (3.33). After integrating over momenta, we have: Setting t =ĉl 2 t /(4π) and recalling the definition of Φ ∞ from (3.82), one trivially derives: Comparing this with the finite volume expression: we can relate the finite and infinite volume expressions for E 4 . In fact, the infinite volume expression can be obtained from the finite volume one by taking theĉ → 0,l → ∞ limit at fixed t , as I fin 8 (t ) vanishes and A(2ĉt ) becomes N 2 −1 N 2 . A detailed discussion on that limit can be found in section 5.
Similar results hold for the other integrals, and thus the infinite volume energy density can be reproduced by performing a simple change in the finite volume formula from eq. (3.46): where I ∞ 9 (t ) = 0 (see appendix C), and I ∞ i (t ) = I i (Φ ∞ , t ) for the rest. Computing the infinite volume integrals in dimensional regularization with d = 4 − 2 , one derives the energy density: which agrees with the result obtained by Lüscher in ref. [3].

't Hooft coupling at one-loop
As we provided in the previous section a regularized expression for the expectation value of the energy density, we are now finally able to focus on several interesting results. Namely, we will in this section derive the running of the coupling, its relation to the MS coupling, obtain the Λ parameter, and present our numerical results for the case of the d t = 2 two-dimensional twist.

Perturbative matching to the MS coupling at one-loop order
Let us begin by recapitulating what has been achieved so far. As we recall, we expanded the observable E(t)/N up to NLO in powers of the 't Hooft coupling: with the leading order term being given by: The function A(x) was defined in eq. (3.41), and the variablesĉ = πc 2 /2 and t = 8t/(cl) 2 were introduced to make the expression more compact.
The NLO contribution is written in terms of twelve integrals given in eqs. (3.62)-(3.73), regulated through analytic continuation in d = 4 − 2 . The leading asymptotic behavior of each integral was identified, and a subtraction procedure was implemented, allowing us to write the energy density at NLO as: Gathering all of these pieces, our results for the expectation value of the energy density can be summarized in the following expression: where C 1 (t) is given by: The perturbative relation to the MS coupling at one-loop order is obtained by simply introducing the expression of the bare coupling in terms of the MS one: leading to: Setting the MS scale to µ = 1/ √ 8t = 1/(cl), the relation at one-loop order between the two couplings reads: λ TGF (l) = λ MS (µ) 1 + c 1 λ MS (µ) , (4.9)

JHEP03(2019)200
with the following matching coefficient at one-loop order: and where we introduced the one-loop constant C 1 : The ratio between Λ parameters in both schemes is then determined, as usual, in terms of the finite one-loop constant c 1 : The purpose of the rest of this section will be to evaluate C 1 numerically, in the case of a single non-trivially twisted plane.

The matching coefficient for a two-dimensional twist
The ingredients required in order to compute the finite constant C 1 , entering the ratio Λ TGF /Λ MS , have been provided in section 3.2. In the specific case of d t = 2, the computational effort that has to be invested in order to determine C 1 is considerably smaller than for d t = 4, as the 8 × 8 matrices entering the expression for Φ are reduced to, at most, 4 × 4. In particular, we have: × Re Θ 0|iB ĉs,ĉu,ĉv,θ − Θ 2 0|iA 0 ĉsl 2 g ,ĉu,ĉvl g , where we defined a 2 × 2 matrix A 0 : as well as a 4 × 4 matrix B containing theθ dependence given by, denoting the twodimensional Levi-Civita symbol: The starting point for the numerical calculation of C 1 will then be given by eqs. (3.89) and (3.93), defining I fin i . All these integrals have been built to be finite, so d can be set to four, and l g to N , in all intervening expressions. The calculation will come in two steps, the first of which will involve using a short Mathematica program to evaluate: for i = 1, · · · 8 and i = 10, · · · 12, (4.17)

JHEP03(2019)200
The required Jacobi theta functions are part of the standard Mathematica package, and for the integration we used the numerical integrators provided by the program by default. The derivatives appearing in some of the integrals were computed using finite differences. The second step is far more complex from a numerical viewpoint, as it encompasses the calculation of: for i = 1, · · · 8 and i = 10, · · · 12, (4.19) Two independent codes were prepared for this second step, one of them written in Mathematica 3 and the other in C++. The former, much like in the first step, made use of the standard Mathematica packages, numerical integrators, and finite differences to compute the integrals, whereas the full details of the inner workings of the latter can be found in appendix D. We will simply mention here that different errors were used for each of the integrals, depending on computation time. The relative errors ranged from 10 −8 in the best cases (for the single integrals), to 10 −3 at worst for I 9 , which was by far the bottleneck.
The value of c also had significant effects, with lower values taking longer times to compute.
Two key aspects are particularly interesting in the analysis of the results: the dependence on c of the coupling at constantθ, and the general dependence inθ.
For an example of the former, we analyzed in detail the case of SU(3) withk = 1, with c ranging from 0.18 to 0.8. The results for C 1 are shown in table 1. Figure 1 displays log(Λ TGF /Λ MS ) as a function of c. In a few points we plot the results obtained with both the Mathematica and the C++ codes, which are perfectly compatible (errors in the data points are smaller than the size of the symbol). The yellow horizontal line shows the result obtained when the gradient flow coupling is evaluated at infinite volume. A detailed analysis on the approach to the infinite volume and the dependence on the number of colors is presented in section 5, but for now we will simply mention that at constant energy scale µ = (cl) −1 and fixed N , taking c → 0 is equivalent to taking the large volume limit, in which log(Λ TGF /Λ MS ) should approach the yellow line in the plot.
As for the study of the general dependence onθ, we considered a series of coprime values ofk and (small) N such thatθ ranged from 0.14 to 0.5. The full results for C 1 are shown in table 2 and figure 2, in which they are plotted as a function ofθ for several values of c. We observe that the dependence onθ is rather smooth for the considered values ofk, N . A discussion about theθ-dependence for larger values of N will be presented in section 5.

Dependence on the number of colors and the magnetic flux
In this section, we will analyze the dependence of λ(cl) on the number of colors N and the angular variableθ =k/l g . We will consider two different limits, both of them taken at fixed value of the renormalized 't Hooft coupling. The first is a singular large N limit in the spirit of those introduced in ref. [39], in which N is sent to infinity while the torus size (8)

JHEP03(2019)200
is sent to zero in such a way as to keepl fixed, and the second is the thermodynamic limit, achieved by simultaneously sending c to zero andl to infinity while keeping cl fixed. The idea that the infinite volume limit can be attained atl → ∞ by sending either the torus size or the number of colors to infinity is implicit in our construction.

Singular large N limit andθ-dependence
Singular large N limits such as the one described above have been employed in various contexts. In ref. [12] the non-perturbative running of the SU(∞) 't Hooft coupling was computed through a step scaling procedure implemented by changing the rank of the gauge group. The calculation was done in the extreme case of TEK reduction on a one-site lattice with an effective size given byl = a √ N , where a denotes the lattice spacing. The continuum limit at fixedl was achieved by sending N to infinity, allowing the authors to compute the evolution of the coupling constant through a wide range of scales, and matching the two-loop perturbative formula at small coupling rather well.
These type of limits have also been considered in the framework of non-commutative field theory. The gauge theory we are considering is equivalent, through the Morita duality, to a non-commutative gauge theory whose rational adimensional non-commutativity parameter is given precisely byθ, a mapping through which the effective torus sizel corresponds directly to the size of the non-commutative torus in the dual theory. One of the proposals raised in ref. [39] was to define non-commutative gauge theories at irrational values ofθ through a sequence of ordinary SU(N i ) twisted Yang-Mills theories with increasing number of colors andθ i =k i /N i →θ. In 2+1 dimensions, ref. [37] has shown that this is only possible, avoiding tachyonic instabilities, for an uncountable zero-measure set of values ofθ, such as for instance a sequence of values ofk and N defined throughk i /N i = F i−2 /F i , where F i denotes the ith term in the Fibonacci sequence. In that case, instabilities in the large N limit are avoided and the limiting sequence tends toθ = (3 − √ 5)/2. In 2+1 dimensions, the condition required to avoid instabilities has been shown to be given in terms of a quantity dubbed Z min : where the symbol ||x|| is used to denote the distance from x to the nearest integer [37,40]. Tachyonic instabilities and symmetry breaking transitions can be avoided as long as Z min > 0.1. Remarkably, this parameter also controls, in 4-dimensional perturbation theory, the size of the contribution of non-planar diagrams to the expectation value of Wilson loops [36]. The limiting procedure to define non-commutative gauge theories at irrational values of the non-commutativity parameter relies on the assumption of continuity inθ. The one-loop matching constant C 1 depends on the choice of the parameter c defining the renormalization scheme, the rank of the group, and the magnetic flux k, and, in particular, given a fixed value of c, one should analyze under which conditions the k and N dependence is fully encoded in the dimensionless ratiok/N definingθ. While a detailed analysis of theθ dependence is beyond the scope of this paper, we did look at the integrals I 1 and I 2 entering the definition of C 1 as representative examples of integrals that are respectively UV divergent and finite after dimensional regularization.      Figures 3 and 4 show how the I 1 and I 2 contributions to C 1 depend onθ for c = 0.15 and c = 0.30 respectively. We have explored many values of N ranging from N = 2 to N = 75025, the latter as part of the aforementioned Fibonacci sequence. For c = 0.3, we noticed that the dependence onθ of both integrals is continuous, with the exception of the point N = 2 in the case of I 1 . As c decreases, however, several other points corresponding to small values of N deviate from the general curve, and, in the case of I 1 , we observe a steep dependence onθ for sequences approaching rational values, in particular fork/N = 0, 1/4, 1/3 and 1/2. A similar dependence onθ has been observed in lattice perturbation theory when considering the contribution at second order of non-planar diagrams to the expectation values of Wilson loops [36], which can be understood in terms of the parameter Z min introduced earlier.
Let us take a look at how the dependence in this Z min quantity enters in the I 1 contribution to C 1 . Theθ-dependent term comes from the function H(s, u, v,θ) defined in eq. (3.78). This contribution is finite in the UV and given by: (n −θ˜ m) 2 + iπ mn .   As all terms included in the sum have a non-zero value ofθ˜ m, UV-finiteness is guaranteed. However, in the limit in which this quantity tends to zero, one would retrieve the divergence present in theθ = 0 term. We will in what follows show that such a limit is approached logarithmically in Z min . Let us begin by considering the leading asymptotic behavior for small x: wheren denotes the integer closest toθ˜ m. Integrating over x, we get: where Z 2 (m) = m 2 ||θm|| 2 . If the argument of the incomplete Γ function is small, this goes as: The logarithmic dependence in Z is tamed by the exponential damping inĉm 2 , but at small enoughĉ this suppression disappears, giving rise to the behavior presented in figure 4a. This is more clearly seen in figure 5 where we show the contribution of I 1 to C 1 as a function of log Z min (N, k). The left plot shows the points for which the minimal value is attained at m = (1, 0), and the right one those with the minimum at m = (2, 0), with the red vertical line in the plots corresponding to Z min = 0.1. Sequences approachingθ = 0 in the left plot andθ = 1/2 in the right one are deep in the region with small Z min , where a tiny change in the value ofθ translates into a large change in the integral. As a final remark, we will point out that the value of Z min stays almost constant along the Fibonacci sequence mentioned earlier, meaning that the results of the integrals will JHEP03(2019)200 depend almost exclusively on the value of c. Therefore, as expected, the singular large N limit can be taken safely along such a sequence, making it optimal, for instance, for the determination of the SU(∞) running coupling using the reduction techniques employed in ref. [12].

Large volume limit
So far, we have been discussing the dependence of the matching constant C 1 on the number of colors and the flux-dependent parameterθ for a fixed value of c, the parameter defining the TGF scheme. In contrast, in this subsection we will be looking at a different type of limit, namely the one in which c tends to zero while the effective size is sent to infinity in such a way as to keep flow time fixed (thus fixing the energy scale of the coupling as well). This limit can be taken in two different ways, either by sending the smallest torus period l to infinity while keeping the rank of the group N fixed, or by sending N to infinity at fixed l. If volume independence holds true, in both cases the infinite volume expression should be recovered, and correspondingly C 1 should vanish. As we recall, at fixed value of t, C 1 is a function of three parameters: c, N and the magnetic flux k. In particular, all of the dependence on the boundary conditions (i.e. the dependence on k) is contained in C 1 , and will vanish in the thermodynamic limit provided C 1 does as well. We will therefore analyze in what follows the behavior of the matching constant in the approach to the thermodynamical limit, along with the size of the finite volume (or finite N ) corrections.
To prepare for such a discussion, we will first take a look at the LO term in the expansion of the energy density, eq. (3.40), with t set to (cl) 2 /8. As we recall, the dependence on c and N came from: where: In the infinite volume limit, understood in the sense of c → 0 at fixed l g , one has F 0 (0, d) = 1 and therefore: leading to a LO term in agreement with the results found in ref. [2]. The leading correction is exponentially suppressed with the square of the volume as: If the large N limit (i.e. large l g ) at fixed l and constant cl g is taken instead, one gets A(πc 2 ) = 1 + O(1/N 2 ), which does indeed correspond to the infinite volume large N limit. The approach to the limit is in that case powerlike, with 1/N 2 corrections. The discussion of the NLO term, on the other hand, is more involved and requires some previous steps to be properly considered. As we recall, the different contributions to JHEP03(2019)200 C 1 can be written in a compact way as: 4 where we used the symbol to denote the integrals appearing in eqs. (3.62)-(3.73) in a generic manner, including the prefactors multiplying the Φ function and derivatives when required. The quantityĤ(s, u, v,θ) is related to the function H(s, u, v,θ) entering the definition of Φ through: 11) and it is given by: (5.12) with: In order to analyze the approach to the infinite volume limit, it is more convenient to look at the expression resulting after Poisson resummation in m. We will, for simplicity's sake, focus on the case of the two-dimensional twist, d t = 2, and will move the full detail of the computations to appendix E for clarity. We will separate each of the contributions to C 1 intoθ−independent andθ−dependent terms, given by: , (5.14) where the functionĤ is obtained by subtracting the zero modes fromĤ after Poisson resummation (see appendix E for the details), and: TD , and depends on two quantities:ĉα andĉαu/s. The simplest case corresponds to integralsĪ 1 ,Ī 2 andĪ 4 , for which bothĉα andĉαu/s tend to zero in theĉ → 0 limit in all of the integration range. The leading contribution, derived in appendix E, is given by: Integrals for which the infinite volume contribution I ∞ i is UV-divergent at d = 4, such asĪ 1 andĪ 4 , have a leading correction that goes as ∼ log(c 2 N 2 ) exp(−1/(cN ) 2 ). I ∞ 2 is UV-finite and the leading correction has a purely exponential decay in the thermodynamic limit, given by exp(−1/(cN ) 2 ). We show in figure 6 the dependence of these integrals on cN for several values ofk and N , plotting their value multiplied by the factor (N 2 − 1) exp(1/(cN ) 2 ). The continuous lines in the plot are given by the formulas presented above and describe very accurately the data for small cN . In the limit obtained by sending N to infinity and c to zero at small, fixed cN , the three integrals also go to zero with corrections of order 1/N 2 .
The general dependence of C 1 on cN as cN → 0 is in fact well described by a formula analogous to eq. (5.20); an example of this for the case of SU(3) is shown in figure 7, where C 1 is displayed as a function of (cN ) 2 . The continuous line in that plot is the result of a fit to the functional form f (cN ) = exp(−1/(cN ) 2 )(α + β log(cN ) + γcN + δc 2 N 2 ). In order to push the calculation of C 1 to smaller values of c, we split it into two pieces, represented by the open blue circles and the yellow squares in the plot. The most relevant part comes from the contributions ofĪ 3 ,Ī 7 ,Ī 10 ,Ī 11 andĪ 12 , which we were able to compute down to values of (cN ) 2 ∼ 0.1. Asymptotically, this piece is described quite well by the function f (cN ), with a leading dependence on c of the form log(cN ) exp(−1/(cN ) 2 ).
In the rest of this section, we will explore how the infinite volume limit is approached for the remaining integrals (excludingĪ 1 ,Ī 2 , andĪ 4 ). The discussion is a bit more complex in their case, as the leading correction goes as c 2 for each of the integrals, but the corrections cancel out when all contributions to C 1 are considered. We will first analyze the case ofĪ 3 in detail to see how the cancellation takes place, and then generalize it to all other cases. For this integral, in theĉ → 0 limit,ĉαu/s goes to zero in the full integration range, and the leading dependence is given by: From this expression one can show (see appendix E for the details) that the dominant correction in the cN → 0 limit is:  with a 1 = −1.76508480122121275 and, for instance, a 2 (N = 3) = 3.59085631503990722. The quantity a 2 (N )/N 2 grows logarithmically with N 2 , as shown in figure 8. One can show that, in the infinite volume limit, all remaining integralsĪ i converge in the same manner, being proportional to I 0 with a proportionality coefficient of +1 for i = 5, 6, 7, of -1 for i = 10, 11, 12 and of 4 and -2 in the cases ofĪ 8 andĪ 9 respectively. Combining eq. (3.46) with these coefficients, it is easy to show that the total contribution of the leading (cN ) 2 term to C 1 vanishes.
We did not analyze in detail how the different integrals approach zero after subtracting the quadratic piece in c, but, based on the results presented in figure 7, we expect other possible power like corrections to cancel out as well when combined to form C 1 , the final result exponentially decaying towards zero with a leading dependence on c of the form ∼ log(cN ) exp −(cN ) −2 /(N 2 − 1). A preliminary analysis was performed for the case ofĪ 3 , with the quantityĪ 3 − I 0 times the factor (N 2 − 1) exp((cN ) −2 ) being shown in figure 9 as a function of cN for several values of N . Each point in that plot was obtained from the exact expression forĪ 3 , and the continuous lines correspond to the approximate expression obtained combining eqs. (5.16) and (5.17). This decomposition is quite useful towards analyzing the N dependence of the integral, and so we displayed each of the two pieces in figures 10a and 10b as a function of the appropriate scaling variable cN .
Theθ-independent term is presented in figure 10a, multiplied by the factor (N 2 − 1) exp((cN ) −2 ) scaling away most of its N dependence. For cN → 0, the integral decays exponentially as ∼ exp(−(cN ) −2 ), whereas in the large N limit at fixed value of cN it goes to zero with quadratic corrections in 1/N 2 . The analysis of theθ-dependent part is more complicated, as one needs to take into account the dependence on the magnetic flux k. The decay ofĪ 3 towards zero is in this case faster than exponential, going as ∼ (cN ) 2 exp(−(cN ) −2 ). This is shown in figure 10b, where we plottedĪ the inverse of this factor times (N 2 − 1) as a function of cN for various values of N and the magnetic flux. In the large N limit taken at fixed cN , this term also scales to zero as 1/N 2 . It would be interesting to study theθ-dependence for large values of N in more detail for both this integral and the others, but such an analysis goes beyond the scope of this paper.

JHEP03(2019)200 6 Summary and conclusions
We computed the perturbative expansion at one-loop order of the SU(N ) twisted gradient flow coupling, including the matching to the MS infinite volume scheme at a renormalization scale µ = 1/(cl) given by a combination of the size of the torus and the rank of the gauge group. The corresponding one-loop finite piece was determined numerically in the case of a two-dimensional non-trivial twist for whichl = lN . The computation was done for a range of values of c (the number relating the energy scale to the size of the torus), of the magnetic flux, and for several values of the rank N of the gauge group, allowing us to obtain the ratio of Λ parameters between the TGF scheme and the MS one.
Moreover, we deemed it interesting to explore the dependence of the coupling on the number of colors and the magnetic flux in a bit more depth, and so we analyzed the dependence of λ TGF in two different limits. First, we studied the limit in which N and the torus size are sent to infinity and zero respectively in such a way as to keepl, and hence the renormalized 't Hooft coupling at scale µ = 1/(cl), fixed. This is a singular large N limit in the spirit of those introduced in [39], albeit a rather non-standard one since non-planar, θ-dependent diagrams survive the limit as long asl is finite. The connection of this case to non-commutative Yang-Mills theory is straightforward through the use of the Morita duality: the non-commutative dual torus is of lengthl and has a dimensionless noncommutativity parameter given byθ =k/N . Our analysis also supports the observation, first presented in [37,40], that the avoidance of tachyonic instabilities when taking the singular limit is only possible for a zero-measure, though uncountable, set of values ofθ. Curiously, one of the successful cases, of limiting parameterθ = (3 − √ 5)/2, relies on a sequence of Fibonacci numbers with k = F i−2 and N = F i with F i denoting the i-th element of the Fibonacci series [40].
The second limit at which we looked was the thermodynamic limit, in which c is sent to zero andl is sent to infinity while keeping the energy scale µ constant. This leads to the one-loop expression of the 't Hooft gradient flow coupling at infinite volume [3]. Our results give support to the reduction idea, in the sense that the SU(∞) coupling in the thermodynamic limit can also be recovered at fixed torus size by sending N , and hencel, to infinity, in which case the limit is approached with 1/N 2 corrections.

A The Feynman rules with twisted boundary conditions
The Feynman rules for the set of irreducible twist tensors used in this work have been derived in various contexts both in the continuum (see for instance [65] and references therein for a review) and in the lattice regularized version of the theory [23,26,33,36]. In this appendix, we will summarize the ones relevant to our work, derived in the continuum.
The set of allowed gauge transformations in our theory will be restricted to those preserving the form of the boundary conditions in eqs. (3.1), (3.2), using the irreducible twist given in eq. (2.11), and the remaining gauge degrees of freedom will be fixed using a generalized covariant gauge of parameter ξ consistent with the boundary conditions. After scaling the gauge potential with the bare coupling g 0 , the Lagrangian density, including the gauge fixing terms, reads: where D µ ≡ ∂ µ + ig 0 A µ is the covariant derivative and c,c denote the ghost fields. One may then obtain the propagators of the gauge and ghost fields using the Fourier expansion of the gauge potential given in eq. (2.14), along with an analogous one for the ghost fields: where the momenta appearing in these expressions are quantized in units of the effective sizel. The Feynman rules for the vertices are then obtained from the commutation relations in eq. (2.17), and are expressed in terms of the momentum-dependent structure constants F (p, q, −q − r). The terms contributing to minus the gauge fixed action are the following: • 3-gluon term: with: • 4-gluon term:

JHEP03(2019)200
with: • Ghost-gluon term: These rules can be easily used to derive different quantities, such as the one-loop correction to the propagator. At order O(g 2 0 ) and in the Feynman gauge (ξ = 1), the vacuum polarization tensor can be obtained as shown in ref. [65]: B Integral form of the energy density at NLO As we recall, the energy density at NLO in the twisted gradient flow scheme can be expressed in terms of several integrals. In section 3.1.3, we chose for both clarity and concision to show a single example of the derivation of these integrals, and left the expression of the full seven O λ 2 0 contributions to the observable E/N in terms of the integrals for this appendix. The contributing E i terms are the following: When summing all of the terms contributing to E (1) (t), the I i terms cancel out, and thus: C Regularization of I 9 We will present here the full details of the procedure to regularize the integral I 9 , defined in eq. (3.70), which differs slightly from the general treatment described in section 3.2.1.
As we recall, the initial integral is split into three terms: with the Heaviside function θ restricting the integration intervals in z. The first term on the r.h.s. of this expression is already finite in four dimensions, whereas the other two will be shown to be so as well after analytical continuation to d = 4. Let us start by discussing the treatment of I 9 (Φ (0) , t ). The original integral can be rewritten as: which, after some algebraic manipulation, becomes: This integral I is in d = 4 − 2 dimensions of order times zero. The asymptotic behavior at z = 0 is then obtained by expanding A(ĉ(2t + z)) around z = 0, leading to: The integral over z presents a pole in 1/ , but it is cancelled when multiplied by I, leading to a final result that is identically zero for d = 4. This I ∞ 9 is precisely the integral appearing in the infinite volume calculation, and vanishes, as we have just seen, for d = 4 in dimensional regularization.
The remaining I 9 (θ(z − 1)Φ (0) , t ) term can be treated in a similar way. One takes the initial expression:

JHEP03(2019)200
and rewrites it in a form identical to eq. (C.2), only with the integral over z restricted to the interval [1, ∞]. After some manipulation, the regularized result becomes: which is finite in d = 4 dimensions, and which we were able to evaluate numerically.

D Numerical implementation of the integration algorithm
In order to perform the numerical computation required to obtain the results presented in section 4.2, we prepared a code in C++ to compute the values of the Φ(s, u, v,θ) functions and their derivatives at any point, and integrate them along the corresponding ranges using the trapezoidal rule up to a target precision. We will in this section begin by explaining how the computation of each Φ(s, u, v) is performed, and then detail how the integration algorithm works.

D.1 Momentum sums
As we recall, we had to compute the following quantity: which was made finite through the procedure explained in section 3.2.1, and where H and Φ (0) were defined as in that very section, taking d = 4. The three r.h.s. terms can be rewritten in terms of momentum sums of the general form: where we introduced a generic matrix X to denote either A 0 or B from eqs. (4.15) and (4.16). We used Poisson resummation to write the sums in terms of both X and its inverse, allowing us to simultaneously compute several equivalent versions of the three terms of Φ fin , which let us exploit the fact that convergence speed depends on the (s, u, v) point being considered to speed up the program. We defined eight quantities to be computed: where we used the shorthands B θ ≡ B(ĉs,ĉu,ĉv,θ), B 0 ≡ B(ĉs,ĉu,ĉv, 0), andB ≡ B(ĉsN 2 ,ĉu,ĉvN, 0) for clarity. Several of these expressions are redundant: allowing us to rewrite the observable Φ fin (s, u, v) in four equivalent forms: We derivated these four equivalent expressions in the integrals in which it was required, simply using the chain rule and computing the derivatives of each E i function when needed. We will skim over the details of the algorithm used to generate the momenta in the sums, simply mentioning that we defined a four-dimensional integer vector M t = (m 1 , n 1 , m 2 , n 2 ) and generated the corresponding combinations of integers m i , n i , using the M → −M symmetry in the integrand to shorten the computation time. The momentum tetrads were generated in an orderly manner, starting with all contributions of the tetrads with |m i |, |n i | = 0, 1, then adding the ones with some |m i |, |n i | = 2, then |m i |, |n i | = 3 and so on and so forth, adding terms with momenta of increasing order until the sum converges (in the sense that we will detail below).
Thus, the code simply runs through momentum tetrads of increasing order, and passes them through a filter that checks whether or not m is proportional to N and whether or not n = 0, computing and adding the relevant exponential terms to each of the eight E i terms. Once every tetrad of a particular order has been processed, the program computes the value of Φ fin up to that order in the four equivalent ways shown earlier, and checks whether the variation of each term between the previous order and the new one is smaller than a set quantity times the value of the function. If that turns out to be the case for any of the four expressions, the sum is considered to have converged and that particular Φ fin is returned as the result. To avoid early spurious convergences, we set a minimum order of four for the sum. The same relative error was also used as the convergence criterion for the integration algorithm (see next subsection), and ranged between 10 −3 and 10 −8 depending on the integral, due to differences in runtime between them.

D.2 Integration algorithm
Now that we explained how the integrand is computed at each point, we may focus on the integration algorithm, for which we chose to use a fairly standard trapezoidal rule for multiple integrals in which the integral along each coordinate is approximated using an increasing number of trapezoids until a target precision is reached. We will begin by quickly illustrating how a generic single-dimensional integral works in our code, generalize it to the multiple ones, and then mention a few specific choices of strategy.
Consider thus a single integral over a finite interval, say for instance the interval z ∈ [0, 1]. The code begins by computing the value of the integrand, which in this case

JHEP03(2019)200
would be the Φ function, at the beginning and end of the interval, and approximates the integral with a trapezoid. The integrand is then determined at the middle point z = 0.5, and the integral is approximated with the two z ∈ [0, 0.5], [0.5, 1] trapezoids. Then, at third order, the integrand is obtained in the midpoints of the previous trapezoids, and the integral is approximated with the four trapezoids [0, 0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 1]. This subdivision generated by computing the integrand at the midpoints goes on, until the variation in the approximated integral between one order and the next is smaller than a set target (the same that we used for the Φ functions above) times the value of the integral at that order, at which point we consider that convergence has been reached and the integral is finished. As we mentioned earlier, in our runs ranged between 10 −3 and 10 −8 .
Multiple integrals are trivial in such a setting: one simply starts with the integral over the outermost coordinate, z, but at every point in which the integrand needs to be determined instead of computing the Φ function, one recursively calls the integration routine to obtain the integral over the next coordinate.
To allow for easier parallelization, and since the integrand tends to have more structure near z = 0, we chose to split the integral in z into a set of pre-chosen subintervals, with a shorter step size at smaller values of z, and treated the integration along each of these subintervals separately. To avoid spurious convergences, we imposed a minimum of eight points in each integration subinterval. Moreover, in the cases in which the integrals went up to infinity in the z coordinate, we ran the integration code up to z max = 10 4 and extrapolated the result by fitting the results of the last ten subintervals to a simple shifted exponential of the form I f it = a 0 − a 1 e −a 2 (z−zmax) , using the fitted a 0 as the final result of the integral. A simple least squares method algorithm was used to perform the fits.
There were a couple of peculiarities worth mentioning regarding integrals I 8 and I 9 . For the former, and after performing a change of variables so that the second integral runs up to x = 1, we noticed that the contribution to the integral is concentrated around z = 0, with the profiles of the integrand over x peaking at small values of z and vanishing after a range ∼ z −1 . This means that the strategy to keep dividing the integration interval into halves in the x coordinate is quite inefficient, as the contribution is concentrated in a small region and one is throwing many points into areas that are effectively zero. To avoid this issue, we chose to subdivide the inner integral into 1,5,50,500 and 1000 equal subintervals as z runs up to 1,10,1000 and 10000 respectively. As soon as the integral over two consecutive subintervals in the x axis vanishes for z > 1, the subintervals that follow are ignored entirely, greatly speeding up the computation without affecting the result.
The case of I 9 is a bit special in that the regularization was different from the other integrals, with a Heaviside θ(1 − z) function being introduced in the integrand (see the end of section 3.2.1 and appendix C for the specifics) and separating the bits before and after z = 1. For the numerical computation, we performed the same change of integral as in I 8 to make the second integral run up to y = 1, but then the Heaviside function became a θ(1 − y z ) function, with the integrands being different before and after this point. As convergence turned out to be painfully slow when both integrands were considered jointly, we simply forced the integrals in y to be split from z = 1 onwards into two subintervals [0, 1/z ] and [1/z , 1], with the convergence of each side being considered separately.

JHEP03(2019)200
Due to the procedure we used to determine the convergence of the integrals, for a given integral I, and dubbing the number of integrals to perform n i (single, double or triple), the final error of the integral is: This comes from the fact that both the error of the Φ functions and the convergence criterion for the integrals is given by the same , so for a single integral: Additional integrals simply add extra 1+ factors, which end up generating the (1+n i ) term. In the cases where the integrals ran up to infinity in z and had to be fitted, we presented as the final error either ∆I or the error from the fit itself, whichever was larger.
Moreover, some issues were caused by some computed quantities hitting machine precision, slowing down the computation while leaving the results effectively unaffected. To deal with them, we introduced several hard cuts in the integrals, integrands and determinants. In particular, we made it so that any Φ function returning a value under 10 −12 , any inner integral returning any value under 5 × 10 −12 (or 10 −10 in the cases of a few intervals in which using 5 × 10 −12 led to severe slowdowns), and any exponential returning a result over 10 −13 is automatically set to be exactly zero. The cut in the integrals is also used in the convergence checks we mentioned earlier: whenever the value of the integral times becomes smaller than the precision cut, the precision cut is used as the convergence criterion instead.
Lastly, we need to mention that, despite the integrals computed being finite, convergence near the point (s, u, v) = (2, 0, 0) can become quite slow, as the integrand approaches machine precision. To address this issue, a cut in u was introduced, setting the integrand to zero when u < 0.01 in the integrals in which such point is part of the integration region (namely, in I i for i = 1, 4,5,7,9,11,12). This cut does not appreciably change the results, as the contribution of the excluded area is well below the uncertainty of the total result. To illustrate this, we show in figure 11 some examples of the profile of the integrand near the aforementioned (s, u, v) = (2, 0, 0) point, in which one can both see that the integrand is indeed finite and that the area excluded by the cut is negligible compared to the rest of the integrand.

E The infinite volume and large N limits
In this appendix, we will derive the formulas mentioned in section 5, which were used to analyze the N andθ dependence of C 1 at NLO in the coupling for the case of a two dimensional twist (d t = 2). As we recall, the contributions to C 1 , barring the one from I 9 which is slightly different, can be written in the form: It is convenient, in order to analyze the infinite volume limit, to look at the expressions resulting after Poisson resummation in m for both theθ-dependent andθ-independent parts. For the latter, Poisson resummation yields: In theθ-dependent case, on the other hand, we begin by rewriting m =ml g + m c , with the components of m c µ taking values in the intervals [−l g /2, l g /2) or [−(l g − 1)/2, (l g − 1)/2] when l g is respectively even or odd. Poisson resummation is then performed with respect tom only, leading to:

JHEP03(2019)200
The termĤ(s, u, v,θ) containing theθ dependence requires a bit more work, but the idea is the same. We begin by rewriting the components of the 4-vector n along the twisted directions as n µ =ñ µ N + n c , with n c a 2-dimensional vector of integers taking values for N even or odd in the respective intervals [−N/2, N/2) or [−(N − 1)/2, (N − 1)/2], and then subtract the terms corresponding to n µ = 0 along periodic directions andñ µ = 0 along the twisted ones. Subtracting the m = 0 terms as well, and adding back once more the doubly subtracted ones, we end up with: iŝ cαu (E.13) where z µ = µν n cν k/N + in cµ v/(ĉN 2 αu).
We may then rewrite each of the integrals contributing to C 1 as the sum of two components I = I TI + I TD , the latter containing all of theθ dependence: and with n c and z µ as defined above. From this expression, one can analyze theĉ → 0 limit, whose approach is driven by two variables:ĉα andĉαu/s. In all contributing integrals but I 8 and I 9 , one of the two variables vanishes for all of the integration range when taking such a limit. The first thing worth noting is the fact that zero modes have already been subtracted from all terms not included in I (0) TI and I (0) TD , and hence the leading order in theĉ → 0 limit for them will be proportional to: 18) which approaches zero at least exponentially in theĉN 2 → 0 limit, and goes, in the large N limit taken keepingĉN 2 constant, as 1/N 2 . In most cases, the leading contribution in theĉ → 0 limit is hence given by I to derive the leading correction to the large volume limit. In the three cases it is given by: 16 3(N 2 − 1) (uα) −2 e − π 2ĉN 2 − e − π cαN 2 − e − πŝ cαuN 2 (E. 19) All three integrals can be analytically approximated with this, leading to: e −(cN ) −2 1 + 3γ E − 3 log 3c 2 N 2 − 3c 2 N 2 , (E.20) e −(cN ) −2 1 − 6c 2 N 2 , (E.21) We will now consider the remaining integrals, looking first at theĉ dependence of I (0) TI and I TD . ForĪ 3 ,Ī 6 andĪ 10 , the variable going to zero in theĉ → 0 limit isĉαu/s, and the leading dependence is given by: whereas forĪ 5 ,Ī 7 ,Ī 11 , andĪ 12 the variable going to zero isĉα, and we have instead: To leading order all these integrals go to zero as ∼ c 2 , with a coefficient depending on N that is identical in absolute value for all of them. We will take a look atĪ 3 as an illustrative example. The leading contribution in thê c → 0 limit for this integral is given by: which allows us to separate A into two parts, one that depends on N and another that does not: Rescaling z to z =ĉz in the first expression and to z =ĉN 2 z in the second, we can decompose the integral into the difference of two piecesĪ 3 , which in the cN → 0 limit become:Ī The leading order result in the cN → 0 limit is thus given by: TD . The cases ofĪ 8 andĪ 9 are shown in the plot as well, which also turn out to be proportional to I 0 with respective coefficients 4 and -2.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.