The twisted gradient flow coupling at one loop

TL;DR

This paper calculates the one-loop running of the SU(N) 't Hooft coupling in a finite volume using twisted boundary conditions, providing a scheme for regularization and scheme matching with potential implications for non-commutative gauge theories.

Contribution

It introduces a method to compute the one-loop running of the coupling in a twisted gradient flow scheme with explicit scheme matching and N dependence analysis.

Findings

Derived the one-loop running of the coupling in a twisted scheme.

Provided a numerical evaluation of the scheme matching coefficient.

Analyzed the N dependence and implications for non-commutative theories.

Abstract

We compute the one-loop running of the $SU(N)$ 't Hooft coupling 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 $\tilde{l}$ 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 $\overline{\rm MS}$ 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 $\Lambda$ 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.

Figures

Tables

Table 1. Table 1 : Results for 𝒞 1 subscript 𝒞 1 \mathcal{C}_{1} for the S ​ U ​ ( 3 ) 𝑆 𝑈 3 SU(3) gauge group and a range of values of c 𝑐 c .

$c$	${𝒞}_1$	$c$	${𝒞}_1$	$c$	${𝒞}_1$	$c$	${𝒞}_1$
0.18	0.224(8)	0.34	0.301(5)	0.50	-1.831(3)	0.66	-4.577(12)
0.19	0.289(7)	0.35	0.228(6)	0.51	-2.014(6)	0.67	-4.727(20)
0.20	0.353(7)	0.36	0.142(2)	0.52	-2.198(3)	0.68	-4.862(17)
0.21	0.404(6)	0.37	0.0530(16)	0.53	-2.383(3)	0.69	-4.998(17)
0.22	0.451(2)	0.38	-0.0464(5)	0.54	-2.569(3)	0.70	-5.119(15)
0.23	0.493(2)	0.39	-0.152(3)	0.55	-2.755(4)	0.71	-5.239(10)
0.24	0.536(3)	0.40	-0.268(2)	0.56	-2.947(8)	0.72	-5.359(10)
0.25	0.557(3)	0.41	-0.385(2)	0.57	-3.125(9)	0.73	-5.460(17)
0.26	0.570(2)	0.42	-0.525(4)	0.58	-3.303(10)	0.74	-5.581(14)
0.27	0.567(4)	0.43	-0.664(3)	0.59	-3.482(10)	0.75	-5.705(10)
0.28	0.558(3)	0.44	-0.813(2)	0.60	-3.646(14)	0.76	-5.806(18)
0.29	0.532(5)	0.45	-0.971(5)	0.61	-3.808(10)	0.77	-5.877(60)
0.30	0.508(4)	0.46	-1.134(5)	0.62	-3.968(12)	0.78	-6.011(40)
0.31	0.473(4)	0.47	-1.302(4)	0.63	-4.125(16)	0.79	-6.141(40)
0.32	0.426(6)	0.48	-1.474(3)	0.64	-4.280(12)	0.80	-6.248(40)
0.33	0.361(3)	0.49	-1.650(3)	0.65	-4.435(16)

Table 2. Table 2 : Results for 𝒞 1 subscript 𝒞 1 \mathcal{C}_{1} for several S ​ U ​ ( N ) 𝑆 𝑈 𝑁 SU(N) gauge groups and values of k ¯ ¯ 𝑘 \bar{k} .

$\overline{k}$	$N$	$\hat{\theta }=\overline{k}/N$	c=0.4	c=0.5	c=0.6	c=0.7	c=0.8
1	7	0.1429	-4.672(15)	-5.814(20)	-6.813(26)	-7.799(39)	-8.693(62)
1	6	0.1667	-4.274(14)	-5.729(19)	-6.979(24)	-8.097(35)	-9.080(47)
1	5	0.2000	-3.417(12)	-5.098(17)	-6.573(23)	-7.811(33)	-8.843(43)
1	4	0.2500	-2.049(12)	-3.808(16)	-5.475(22)	-6.833(30)	-7.912(40)
2	7	0.2857	-1.187(13)	-2.891(15)	-4.634(20)	-6.050(29)	-7.156(39)
1	3	0.3333	-0.261(14)	-1.818(14)	-3.614(19)	-5.087(29)	-6.220(38)
3	8	0.3750	0.327(14)	-1.073(16)	-2.888(19)	-4.395(26)	-5.545(37)
2	5	0.4000	0.583(14)	-0.724(16)	-2.542(19)	-4.064(26)	-5.222(37)
3	7	0.4286	0.791(12)	-0.418(16)	-2.236(19)	-3.771(26)	-4.937(36)
5	11	0.4545	0.911(09)	-0.228(15)	-2.045(19)	-3.587(26)	-4.757(36)
1	2	0.5000	1.077(13)	-0.092(16)	-1.914(19)	-3.461(26)	-4.634(36)

Equations

$${\partial }_t{B}_{\mu }\left(x,t\right)=D_{\nu }{G}_{\nu \mu }\left(x,t\right),B_{\mu }\left(x,0\right)=A_{\mu }\left(x\right),$$

$$D_{\mu }{B}_{\nu }\left(x,t\right)$$

$$G_{\mu \nu }\left(x,t\right)$$

$$\langle E(t)\rangle=\frac{1}{2}\Big{\langle}\mathrm{Tr}\Big{(}G_{\mu\nu}(x,t)G_{\mu\nu}(x,t)\Big{)}\Big{\rangle}.$$

$$\lambda (\mu )={\frac{128{\pi }^2{N}^2}{(d-1)(N^2-1)}⟨\frac{t^{2}E(t)}{N}⟩|}_{t=1/(8{\mu }^2)},$$

$${\lambda }_{\mathrm{TGF}}(c\text{tl})={\mathcal{F}(c)⟨\frac{t^{2}E(t)}{N}⟩|}_{t=c^2{\text{tl}}^2/8},$$

$$A_{\mu }(x+l_{\nu }\hat{\nu })={\Gamma }_{\nu }{A}_{\mu }(x){\Gamma }_{\nu }^{†},$$

$${\Gamma }_{\mu }{\Gamma }_{\nu }=Z_{\mu \nu }{\Gamma }_{\nu }{\Gamma }_{\mu },$$

$$Z_{\mu \nu }=\exp \{i2\pi n_{\mu \nu }/N\}.$$

$$n_{\mu \nu }={ϵ}_{\mu \nu }\frac{kN}{l_g},$$

$$Z_{\mu \nu }=\exp \left\{i2\pi {ϵ}_{\mu \nu }\frac{k}{l_g}\right\}.$$

$$\hat{\Gamma }(q)=\frac{1}{\sqrt{2N}}{e}^{i\beta \left(q\right)}{\Gamma }_0^{s_0\left(q\right)}\dots {\Gamma }_{d_t-1}^{s_{d_t-1}\left(q\right)},$$

$${\Gamma }_{\nu }\hat{\Gamma }(q){\Gamma }_{\nu }^{†}=e^{i{q}_{\nu }{l}_{\nu }}\hat{\Gamma }(q),$$

$$A_{\mu }(x)=V^{-\frac{1}{2}}\underset{q}{\sum ^{\prime }}{e}^{iqx}{\hat{A}}_{\mu }(q)\hat{\Gamma }(q),$$

$$s_{\mu }(q)={\widetilde{ϵ}}_{\mu \nu }\bar{k}\frac{l_g{l}_{\nu }{q}_{\nu }}{2\pi }\in \mathbb{Z},$$

$$k\bar{k}=1(\mathrm{mod}{l}_g),\sum _{\nu }{\widetilde{ϵ}}_{\mu \nu }{ϵ}_{\nu \rho }={\delta }_{\mu \rho }.$$

$$\left[\hat{\Gamma }\left(p\right),\hat{\Gamma }\left(q\right)\right]=iF\left(p,q,-p-q\right)\hat{\Gamma }\left(p+q\right),$$

$$F\left(p,q,-p-q\right)=-\sqrt{\frac{2}{N}}\sin \left(\frac{1}{2}{\theta }_{\mu \nu }{p}_{\mu }{q}_{\nu }\right),$$

$${\theta }_{\mu \nu }=\frac{{\text{tl}}_{\mu }{\text{tl}}_{\nu }}{2\pi }{\widetilde{ϵ}}_{\mu \nu }\hat{\theta },\hat{\theta }=\frac{\bar{k}}{l_g}.$$

$$A_{\mu }(x+l\hat{\nu })$$

$$A_{\mu }(x+\text{tl}\hat{\nu })$$

$$\Omega (x+l\hat{\nu })$$

$$\Omega (x+\text{tl}\hat{\nu })$$

$${\partial }_t{B}_{\mu }\left(x,t\right)=D_{\nu }{G}_{\nu \mu }\left(x,t\right)+\xi D_{\mu }{\partial }_{\nu }{B}_{\nu }\left(x,t\right),B_{\mu }\left(x,0\right)=g_0{A}_{\mu }\left(x\right),$$

$$B_{\mu }(x,t)=\sum _k{g}_0^k{B}_{\mu }^{(k)}(x,t),B_{\mu }^{(k)}(x,0)={\delta }_{k1}{A}_{\mu }(x).$$

$$B_{\mu }^{(k)}(x,t)=V^{-\frac{1}{2}}\underset{q}{\sum ^{\prime }}{e}^{iqx}{\hat{B}}_{\mu }^{(k)}(q,t)\hat{\Gamma }(q).$$

$$E(t)$$

$$+2i{g}_0^3\mathrm{Tr}\left({\partial }_{\mu }{B}_{\nu }^{(1)}\left[B_{\mu }^{(1)},B_{\nu }^{(1)}\right]\right)$$

$$+2g_0^3\mathrm{Tr}\left({\partial }_{\mu }{B}_{\nu }^{(1)}{\partial }_{\mu }{B}_{\nu }^{(2)}-{\partial }_{\nu }{B}_{\mu }^{(1)}{\partial }_{\mu }{B}_{\nu }^{(2)}\right)$$

$$+g_0^4\mathrm{Tr}\left({\partial }_{\mu }{B}_{\nu }^{(2)}{\partial }_{\mu }{B}_{\nu }^{(2)}-{\partial }_{\mu }{B}_{\nu }^{(2)}{\partial }_{\nu }{B}_{\mu }^{(2)}\right)$$

$$-\frac{1}{2}{g}_0^4\mathrm{Tr}\left({\left[B_{\mu }^{(1)},B_{\nu }^{(1)}\right]}^2\right)$$

$$+2i{g}_0^4\mathrm{Tr}\left({\partial }_{\mu }{B}_{\nu }^{\left(2\right)}\left[B_{\mu }^{(1)},B_{\nu }^{(1)}\right]+{\partial }_{\mu }{B}_{\nu }^{(1)}\left[B_{\mu }^{(1)},B_{\nu }^{(2)}\right]+{\partial }_{\mu }{B}_{\nu }^{(1)}\left[B_{\mu }^{(2)},B_{\nu }^{(1)}\right]\right)$$

$$+2g_0^4\mathrm{Tr}\left({\partial }_{\mu }{B}_{\nu }^{(3)}{\partial }_{\mu }{B}_{\nu }^{(1)}-{\partial }_{\mu }{B}_{\nu }^{(3)}{\partial }_{\nu }{B}_{\mu }^{(1)}\right)+\mathcal{O}(g_0^5).$$

$${\mathcal{E}}_0(t)$$

$${\mathcal{E}}_1(t)$$

Full text

††preprint:

IFT-UAM/CSIC-19-26

1 Introduction

In recent years, the continuous smoothing procedure known as the gradient flow  Narayanan:2006rf ; Luscher:2009eq ; Luscher:2010iy 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 $\Lambda$ parameter Brida:2014joa ; DallaBrida:2016kgh ; Bruno:2017gxd ; Korzec:2017ypb ; Ishikawa:2017xam to the study of Yang-Mills theories with near conformal behavior Lin:2015zpa ; Fodor:2016zil ; Hasenfratz:2016dou , or with large number of colors GarciaPerez2015 . Several coupling renormalization schemes based on gradient flow techniques have been proposed to that end Luscher:2010iy ; Fodor:2012td ; Borsanyi2012 ; Fritzsch:2013je ; Ramos:2014kla . 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) Luscher:2010iy and next-to-next-to-leading order (NNLO) Harlander:2016vzb , while for finite volume it has been done up to NNLO in the Schrödinger functional scheme using numerical stochastic perturbation theory DallaBrida:2016dai ; DallaBrida:2017tru .

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. Ramos:2014kla . 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) tHooft:1979rtg . 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_{\mu}=0$ present with PBC GonzalezArroyo:1981vw . The usefulness of TBC for perturbation theory was first formulated in the context of volume reduction in large $N$ Yang-Mills theory GonzalezArroyo:1982ub ; GonzalezArroyo:1982hz , and was then extended to various other contexts at finite and large $N$ Fabricius:1985jw ; Luscher:1985zq ; Luscher:1985wf ; Coste:1986cb ; Hansson:1986ia ; GonzalezArroyo:1988dz ; Daniel:1989kj ; Daniel:1990iz ; Snippe:1996bk ; Snippe:1997ru ; Perez:2013dra ; Gonzalez-Arroyo:2014dua ; Perez:2017jyq ; Perez:2018afi ; Gonzalez-Arroyo:2019zfm . Despite these advantanges, 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 GarciaPerez2015 , and set the $SU(N)$ running coupling scale to be proportional to an effective size $\tl$ 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 $\mu=1/(c\tl)$, we will consider two different types of $N\rightarrow\infty$ limits, a thermodynamic limit in which the effective size is sent to infinity as $c$ goes to zero while $\mu$ is kept constant, and a second one, called singular in AlvarezGaume:2001tv , in which $N$ is sent to infinity as the torus size is shrunk to zero at constant $c$ in such a way as to keep $\tl$ fixed Perez:2013dra ; GarciaPerez2015 ; Chamizo:2016msz ; Perez:2018afi  111This type of limit has also been considered in other works such as Guralnik:2001pv ; Unsal:2008ch ; Unsal:2010qh . See Unsal:2018 and references therein as well for more examples.. In ref. GarciaPerez2015 this last limit was used to compute the $SU(\infty)$ 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 Eguchi:1982nm , 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 Bhanot:1982sh , several alternative proposals were formulated Bhanot:1982sh ; Gross:1982at ; Narayanan:2003fc ; Kovtun:2007py ; Unsal:2008ch , one of which was the use of twisted boundary conditions GonzalezArroyo:1982hz ; GonzalezArroyo:1982ub , which has proven very successful provided the twist tensor is judiciously chosen GonzalezArroyo:2010ss ; GonzalezArroyo:2012fx ; Gonzalez-Arroyo:2013bta ; Gonzalez-Arroyo:2014dua ; Perez:2015yna ; Gonzalez-Arroyo:2015bya ; Perez:2017jyq ; Gonzalez-Arroyo:2018aus . The idea of volume reduction with TBC was extended to the continuum theory in ref. GonzalezArroyo:1983ac , constituting the first formulation of the Feynman rules of a non-commutative Yang-Mills theory Douglas:2001ba . The singular limit in ref. AlvarezGaume:2001tv was first formulated precisely in the context of such non-commutative theories, with the effective torus size $\tl$ 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 Chamizo:2016msz ; Perez:2018afi 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 Guralnik:2002ru , 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 sec. 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 size $\tl$), along with some specifics about the implementation of TBC in our setting. In particular, we detailed the orthogonal twist used throughout the paper. Sec. 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 Luscher:2010iy , 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 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 sec. 4, which contains both the matching to the $\MS$ scheme and the ratio of $\Lambda$ 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. Sec. 5 discusses the dependence of the coupling on the number of colors, following similar arguments to those in Perez:2018afi ; Perez:2017jyq . A summary of results is presented in sec. 6. Many technicalities were moved for clarity to appendices A- D, including details on the algorithms used to compute the $\Lambda$ parameter.

2 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$, $\langle t^{2}E(t)\rangle$ is a renormalized quantity and, at leading order in perturbation theory, is proportional to the $\MS$ coupling at a scale $\mu=1/\sqrt{8t}$, which leads to a natural definition of a renormalized coupling constant Luscher:2010iy . 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. Ramos:2014kla , but differs slightly from it for reasons that will become clear in what follows.

2.1 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_{\mu}\left(x,t\right)$ matching the Yang-Mills one $A_{\mu}\left(x\right)$ 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:

[TABLE]

where $D_{\mu}$ and $G_{\mu\nu}$ respectively stand for the covariant derivative and field strength tensor of the flow fields:

[TABLE]

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 Luscher:2011bx .

The renormalized gradient flow coupling can then be defined in terms of the energy density of the flowed field:

[TABLE]

In infinite volume, this quantity can be used to define a renormalized $SU(N)$ ’t Hooft coupling at an energy scale $\mu$, given by Luscher:2010iy :

[TABLE]

where $d$ stands for the number of space-time dimensions of the theory. Finite volume gradient flow schemes Fodor:2012td ; Borsanyi2012 ; Fritzsch:2013je ; Ramos:2014kla ; GarciaPerez2015 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 $\mu=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 in the finite-volume schemes proposed in refs. Ramos:2014kla ; GarciaPerez2015 ; Keegan:2015lva . 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 tHooft:1979rtg , and such that the torus has a period $l$ in $d_{t}$ (twisted) directions, and $\tl=lN^{2/d_{t}}$ 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 GarciaPerez2015 :

[TABLE]

where $\mathcal{F}(c)$ is a constant defined in such a way as to have $\lambda_{\rm TGF}(\tilde{l})=\lambda_{0}+{\mathcal{O}}(\lambda_{0}^{2})$, in terms of the bare ’t Hooft coupling $\lambda_{0}$.

2.2 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 GarciaPerez2015 ; Keegan:2015lva . 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 of $\tl$ as the relevant length scale.

Let us start by considering a $SU(N)$ gauge theory defined on a $d$-dimensional torus of length $l_{\mu}$ in each direction, and focus our analysis in the specific case of four dimensions for a gauge potential that satisfies ’t Hooft TBC tHooft:1979rtg . 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:

[TABLE]

where $\Gamma_{\nu}$ are four $SU(N)$ matrices known as twist eaters which satisfy:

[TABLE]

with $Z_{\mu\nu}$ an element of the center of the gauge group, written in terms of an antisymmetric tensor of integers $n_{\mu\nu}$ as:

[TABLE]

The twist tensor $n_{\mu\nu}$ is preserved under gauge transformations, and uniquely characterizes the boundary conditions. It is said to be orthogonal when $\kappa(n)=\epsilon_{\mu\nu\rho\sigma}n_{\mu\nu}n_{\rho\sigma}/8=0$ (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 $\Gamma_{\mu}$ 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 tHooft:1981nnx ; vanBaal:1983eq 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 ga:torus . 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. (8) is satisfied. We will write our orthogonal twist tensor in the form:

[TABLE]

where $l_{g}=N^{2/d_{t}}$ 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 $\epsilon_{01}=-\epsilon_{10}=1$, and by choosing $\epsilon_{\mu\nu}=0$ in any other plane, whereas to twist all planes non-trivially one must instead take $d_{t}=4$, and set $\epsilon_{\mu\nu}$ to be antisymmetric and equal to 1 whenever $\mu<\nu$. With this choice:

[TABLE]

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 GonzalezArroyo:1982hz , as one can see for instance in Perez:2017jyq for the general treatment when the torus is discretized on a lattice, or in Perez:2013dra 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:

[TABLE]

with $s_{\mu}(q)\in\mathbb{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_{\mu}(q)=0\,({\rm 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 $\hat{\Gamma}(q)$ satisfies:

[TABLE]

with no summation over $\nu$ implied, the boundary conditions in (7) are trivially implemented through the Fourier expansion:

[TABLE]

where $V\equiv\prod_{\mu}l_{\mu}$, and the prime in the sum denotes the exclusion of the momenta for which $\hat{\Gamma}(q)\propto\mathbb{I}$. In the periodic directions, for which $\Gamma_{\nu}\propto\mathbb{I}$, the momenta are as usual quantized in units of $2\pi/l_{\nu}$. This is however not the case for the twisted directions, where a solution is provided by:

[TABLE]

where $\bar{k}$ and $\tilde{\epsilon}_{\mu\nu}$ are given by:

[TABLE]

The momentum along the twisted $\mu$ directions is thus quantized in units of $2\pi/\tl_{\mu}$, with $\tl_{\mu}\equiv l_{\mu}l_{g}$. For this choice of $s_{\mu}$, the group structure constants in the $\hat{\Gamma}(q)$ basis become momentum dependent and are given by:

[TABLE]

with

[TABLE]

and

[TABLE]

The tracelessness of the $\hat{\Gamma}(q)$ matrices thus forbids momenta such that $\tl_{\mu}q_{\mu}=0\,({\rm mod}\,2\pi 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, $\tl_{\mu}=l_{\mu}l_{g}=l_{\mu}N^{2/d_{t}}$, for the latter. This observation has led us to a specific choice of torus size to define the TGF coupling in eq. (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 and $\tl$ 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\pi/\tl$, and we will use this effective size $\tl$ as the renormalization scale for the running coupling.

3 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. Luscher:2010iy . 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 Brezin:1985xx – in a way that will be detailed in this section.

3.1 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:

[TABLE]

where $\Gamma_{\mu}$ satisfies eq. (8), and with a twist tensor of the form shown in eq. (11). This restricts the set of allowed gauge transformations $\Omega(x)$ down to those preserving the form of the twist matrices, i.e. those satisfying:

[TABLE]

These boundary conditions are implemented through the Fourier expansion of the gauge field given in eq. (14). In the specific case of the asymmetric torus that we are considering, the torus volume is given by $V=l^{d_{t}}\tl^{\,d-d_{t}}$, and momenta in all directions are quantized in terms of the effective size $\tl$. As we recall, the prime in the sum in eq. (14) denotes the exclusion of all momenta for which $\tl q_{\mu}=0$ (mod $2\pi 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_{\mu}=0$, zero-action solution. We start in $d=4-2\epsilon$ dimensions by scaling the original gauge potential with the bare coupling, $A_{\mu}(x)\rightarrow g_{0}A_{\mu}(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:

[TABLE]

$\xi$ 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 Luscher:2010iy , 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 (20) at any given flow time Ramos:2014kla . These modified flow equations can be solved order by order in $g_{0}$ by expanding the flow field in powers of the coupling:

[TABLE]

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:

[TABLE]

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_{\mu}^{(n)}\equiv B_{\mu}^{(n)}(x,t)$, one gets, up to order $g_{0}^{4}$:

[TABLE]

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 $\hat{\Gamma}(q)$ matrices – see eqs. (17), (18). For the sake of completeness we give below the seven different terms contributing to the expectation value of $\langle E(t)\rangle$ arising at order $g_{0}^{4}$, with an additional $1/N$ normalization factor added for later convenience. Each term can be identified with one of the lines in eq. (27):

[TABLE]

[TABLE]

The shorthand notation $\sum p_{i}$ in the $\delta$ functions was used to denote the sum over all present momenta for each term. The ${\mathcal{E}}_{0}$ term will turn out to be a combination of a leading $\mathcal{O}\left(g_{0}^{2}\right)$ term and an $\mathcal{O}\left(g_{0}^{4}\right)$ correction, whereas all other terms will turn out to be $\mathcal{O}\left(g_{0}^{4}\right)$.

Then, the next step is to relate the flow fields to the actual gauge fields $A_{\mu}\left(x\right)$, for which we will need to obtain an order-by-order solution to the flow equations.

3.1.1 Solving the flow equations in the TGF scheme

Let us consider the flow equation (24) with the gauge parameter $\xi$ set to unity. This was already solved by Lüscher for the infinite volume case Luscher:2010iy , but the results in finite volume are slightly different. Expanding the fields in perturbation theory, as in eq. (25), 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:

[TABLE]

The first three orders will be enough to obtain the observable at order $\mathcal{O}\left(g_{0}^{4}\right)$:

[TABLE]

We may define a momentum space version of $R_{\mu}^{\left(i\right)}$:

[TABLE]

under which the $R_{\mu}$ terms read:

[TABLE]

[TABLE]

In terms of $\hat{B}_{\mu}\left(q,t\right)$, the flow equation in momentum space becomes:

[TABLE]

whose solution is immediate at first order:

[TABLE]

and which can be solved for the next two orders by directly integrating $R_{\mu}^{\left(i\right)}$:

[TABLE]

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. (28)- (34) in terms of sums over momenta. Introducing for the sake of readability the symbol:

[TABLE]

and after quite a bit of algebra, we ended up with:

[TABLE]

[TABLE]

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 ${\cal O}(\lambda_{0}^{2})$ terms have been discarded, as the $F^{2}(q,r,-q-r)$ factors automatically vanish for such momenta.

3.1.2 The energy density at LO

As we recall, our aim in this paper was to obtain a perturbative expansion of the observable $\langle E(t)/N\rangle$ at NLO, which in powers of the bare ’t Hooft coupling can be parametrized as:

[TABLE]

We will begin by deriving the leading order term from the formulas in the previous subsection, given by:

[TABLE]

It will be convenient to introduce a few auxiliary variables and functions:

[TABLE]

in terms of which we may write:

[TABLE]

This ${\mathcal{A}}$ function can be expressed in terms of Jacobi theta functions $\theta_{3}$ as:

[TABLE]

where we used Poisson resummation to rewrite the theta functions:

[TABLE]

The leading order infinite volume expression is retrieved in the $c\rightarrow 0$, $\tl\rightarrow\infty$ limit, taken in such a way as to keep $c\tl$ fixed. In that limit:

[TABLE]

leading to:

[TABLE]

in agreement with the results in ref. Luscher:2010iy .

3.1.3 The energy density at NLO

As for the subleading ${\cal O}(\lambda_{0}^{2})$ term coming from eqs. (48) - (54), 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:

[TABLE]

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, ${\mathcal{E}}_{4}$ in eq. (52), and show the remaining $\mathcal{E}_{i}$ contributions in terms of the basic integrals in appendix B. The ${\mathcal{E}}_{4}$ contribution is given by:

[TABLE]

Using Schwinger parametrization to lift the momenta from the denominator, we then defined an integral:

[TABLE]

so as to rewrite:

[TABLE]

with the structure constants entering this expression through the definition of the symbol $\underset{q,r}{\widehat{\sum}}$ given in eq. (46).

The presence of the structure constants in each ${\mathcal{E}}_{i}$ will let us formulate the integrands in terms of Siegel theta functions. We may indeed rewrite $NF^{2}$ as:

[TABLE]

a substitution under which a generic integrand of the form:

[TABLE]

becomes:

[TABLE]

where we rescaled the variables $s\equiv 4\pi\tl^{-2}\tilde{s}$, $u\equiv 4\pi\tl^{-2}\tilde{u}$, $v\equiv 4\pi\tl^{-2}\tilde{v}$, and where we used the quantization of momenta in the twisted finite box to rewrite $q$ and $r$ in terms of integers.

The connection to Siegel theta functions becomes clear by introducing the function

[TABLE]

with

[TABLE]

In this expression $\mathbb{I}_{d}$ denotes the $d\times d$ identity matrix and the sum over $M$ denotes the sum over the corresponding integers $m_{\mu},n_{\nu}$, regrouped into a $2d$-dimensional column vector. Recalling the definition of the Siegel theta functions:

[TABLE]

this matricial expression takes the form:

[TABLE]

Using this notation, the integral entering ${\mathcal{E}}_{4}$ reads:

[TABLE]

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^{\prime}$ and ${\hat{c}}$ defined in eq. (57) and rescaling $z$ appropriately we have:

[TABLE]

Introducing an auxiliary function $\Phi(s,u,v,\hat{\theta})$ that incorporates the normalization factor in front of the integral:

[TABLE]

with:

[TABLE]

we rewrote the integrals in a fairly basic form allowing us to evaluate them numerically. For instance, for the integral in ${\mathcal{E}}_{4}$:

[TABLE]

A similar procedure can be followed for all the terms contributing to the energy density at NLO, leading to the result in eq. (65) for ${\mathcal{E}}^{(1)}(t)$, where the twelve intervening integrals are:

[TABLE]

3.2 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. Luscher:2010iy 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 $\Phi$ function:

[TABLE]

The real part of the matrix $A({\hat{c}}s,{\hat{c}}u,{\hat{c}}v,\hat{\theta})$, obtained by setting $\hat{\theta}=0$ in eq. (77), is a positive definite symmetric matrix as long as $\det A({\hat{c}}s,{\hat{c}}u,{\hat{c}}v,\hat{\theta}=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:

[TABLE]

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$. 222In some cases divergences happen for the points $(s,u,v)=(2t^{\prime},2t^{\prime},2t^{\prime})$, but, using a particular momentum shift, they can be moved to $(s,u,v)=(2t^{\prime},0,0)$. Indeed, the substitution: $u^{\prime}=s+u-2v,\,v^{\prime}=v-s$ can be implemented by shifting $m\rightarrow m-n$ within the momentum sums. Such a shift had already been applied in order to give the definitions of the integrals in the previous subsection.

The analysis of the asymptotic behavior of the integrals is much clearer once we apply Poisson resummation as in eq. (62) to each component of $n$ in the definition of $\Theta$:

[TABLE]

Whenever $\hat{\theta}\tilde{\epsilon}m\notin\mathbb{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 $\hat{\theta}$ 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\rightarrow 0$, as:

[TABLE]

This observation allows us to isolate the asymptotic divergence by identifying the cases for which $\hat{\theta}\tilde{\epsilon}m\in\mathbb{Z}^{d}$. The first case in which this occurs is whenever $\hat{\theta}\equiv\bar{k}/l_{g}=0$, for any value of $m$. For nonzero $\hat{\theta}$, it will happen whenever $\tilde{\epsilon}m=0$ (mod $l_{g})$. Since the vector $\tilde{\epsilon}m$ has nonvanishing components only along the twisted directions, this will be the case whenever $m_{\mu}=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 $\hat{\theta}=0$, terms with $n=0$ and any value of $m$.

• •

For $\hat{\theta}\neq 0$, terms with $n=0$ (after shifting away the $\hat{\theta}$ dependence), and $m_{\mu}=0\text{ (mod }l_{g})$ in all twisted directions at once.

With this, we may begin the discussion on how the divergent integrals can be regularized.

3.2.1 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,\hat{\theta})$ given by:

[TABLE]

with the usual meaning for the prime in the sum over $m$. The $\Phi$ function entering the integrals can then be rewritten as:

[TABLE]

which is quite advantageous, as the explicit exclusion of the momenta $m$ proportional to $l_{g}$ from the sum automatically makes the term in $\hat{\theta}\neq 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:

[TABLE]

with $\alpha$ defined as in eq. (98), 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 $\Phi^{(0)}$ in terms of the function ${\mathcal{A}}(x)$ from eq. (58):

[TABLE]

For reasons that will become clear later, we define:

[TABLE]

in terms of which we may rewrite:

[TABLE]

This formulation will be useful to analyze the asymptotic UV behavior of the integrals resulting from replacing the original function $\Phi$ in the integrand by $\Phi^{(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:

[TABLE]

which, substituting in the expression for $\Phi^{\infty}$, yields:

[TABLE]

The asymptotic behavior at small $x$ can then be obtained by expanding ${{\mathcal{A}}}\left(2{\hat{c}}t^{\prime}-{\hat{c}}x/2\right)$ 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:

[TABLE]

Notice that in this expression the entire momentum dependence has been factorized into the normalization constant ${{\mathcal{A}}}(2{\hat{c}}t^{\prime})$, which happens to be the same factor that appeared at leading order – see eq. (59). The integral $I_{1}(\Phi^{\infty},t^{\prime})$ can then be evaluated in dimensional regularization with $d=4-2\epsilon$, leading to:

[TABLE]

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 ${\mathcal{A}}({\hat{c}}\alpha)$ appearing in the definition of $\Phi^{(0)}$ around $u=0$, retain the leading term, and then use it to define:

[TABLE]

Remarkably, the integrals $I_{i}(\Phi^{\infty},t^{\prime})$ match the ones appearing in the infinite volume calculation (up to a factor depending on $N$), which we will present in sec. 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:

[TABLE]

and another one, shown in eq. (111) 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/\epsilon$ 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:

[TABLE]

with the Heaviside step function $\theta$ 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_{9}^{\text{reg}}$ these analytic continuations, we end up with:

[TABLE]

And therefore:

[TABLE]

and:

[TABLE]

3.3 Infinite volume limit

The expression of the energy density in infinite volume can be easily retrieved (see Perez:2017jyq ) by making the following substitutions in eqs. (47) - (54):

[TABLE]

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, much like what happened in the finite volume case. We will first present the case of ${\mathcal{E}}_{4}$ as an illustrative example, and then present the results for the general case.

The infinite volume expression for ${\mathcal{E}}_{4}$ is obtained making the substitutions from eqs. (118) and  (119) in eq. (52). After integrating over momenta, we have:

[TABLE]

Setting $t={\hat{c}}\tl^{2}t^{\prime}/(4\pi)$ and recalling the definition of $\Phi^{\infty}$ from (105), one trivially derives:

[TABLE]

Comparing this with the finite volume expression:

[TABLE]

we can relate the finite and infinite volume expressions for ${\mathcal{E}}_{4}$. In fact, the infinite volume expression can be obtained from the finite volume one by taking the $\hat{c}\rightarrow 0$, $\tilde{l}\rightarrow\infty$ limit at fixed $t^{\prime}$, as $I_{8}^{\text{fin}}(t^{\prime})$ vanishes and $\mathcal{A}(2{\hat{c}}t^{\prime})$ becomes $\frac{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. (65):

[TABLE]

where $I_{9}^{\infty}(t^{\prime})=0$ (see appendix C), and $I_{i}^{\infty}(t^{\prime})=I_{i}(\Phi^{\infty},t^{\prime})$ for the rest.

Computing the infinite volume integrals in dimensional regularization with $d=4-2\epsilon$, one derives the energy density:

[TABLE]

which agrees with the result obtained by Lüscher in ref. Luscher:2010iy .

4 ’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 $\Lambda$ parameter, and present our numerical results for the case of the $d_{t}=2$ two-dimensional twist.

4.1 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 $\langle E(t)/N\rangle$ up to NLO in powers of the ’t Hooft coupling:

[TABLE]

with the leading order term being given by:

[TABLE]

The function ${\mathcal{A}}(x)$ was defined in eq. (60), and the variables ${\hat{c}}=\pi c^{2}/2$ and $t^{\prime}=8t/(c\tl)^{2}$ were introduced to make the expression more compact.

The NLO contribution is written in terms of twelve integrals given in eqs. (85)- (96), regulated through analytic continuation in $d=4-2\epsilon$. 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:

[TABLE]

All of the $1/\epsilon$ poles arising in dimensional regularization are contained in ${\mathcal{E}}_{\rm div}^{(1)}(t)$, a quantity that can be trivially rewritten in terms of the infinite volume result ${\mathcal{E}}_{\infty}^{(1)}(t)$ as:

[TABLE]

Gathering all of these pieces, our results for the expectation value of the energy density can be summarized in the following expression:

[TABLE]

where $C_{1}(t)$ is given by:

[TABLE]

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:

[TABLE]

leading to:

[TABLE]

Setting the $\MS$ scale to $\mu=1/\sqrt{8t}=1/(c\tl)$, the relation at one-loop order between the two couplings reads:

[TABLE]

with the following matching coefficient at one-loop order:

[TABLE]

and where we introduced the one-loop constant $\mathcal{C}_{1}$:

[TABLE]

The ratio between $\Lambda$ parameters in both schemes is then determined, as usual, in terms of the finite one-loop constant $c_{1}$:

[TABLE]

The purpose of the rest of this section will be to evaluate $\mathcal{C}_{1}$ numerically, in the case of a single non-trivially twisted plane.

4.2 The matching coefficient for a two-dimensional twist

The ingredients required in order to compute the finite constant $\mathcal{C}_{1}$, entering the ratio $\Lambda_{\rm TGF}/\Lambda_{\MS}$, have been provided in sec. 3.2. In the specific case of $d_{t}=2$, the computational effort that has to be invested in order to determine $\mathcal{C}_{1}$ is considerably smaller than for $d_{t}=4$, as the $8\times 8$ matrices entering the expression for $\Phi$ are reduced to, at most, $4\times 4$. In particular, we have:

[TABLE]

where we defined a $2\times 2$ matrix $A_{0}$:

[TABLE]

as well as a $4\times 4$ matrix $B$ containing the $\hat{\theta}$ dependence given by, denoting $\epsilon$ the two-dimensional Levi-Civita symbol:

[TABLE]

The starting point for the numerical calculation of $\mathcal{C}_{1}$ will then be given by eqs. (112) and  (116), defining $I_{i}^{\rm fin}$. 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:

[TABLE]

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:

[TABLE]

Two independent codes were prepared for this second step, one of them written in Mathematica 333In this case it turned out to be convenient to evaluate the two steps of the calculation jointly. 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 $\hat{\theta}$, and the general dependence in $\hat{\theta}$.

For an example of the former, we analyzed in detail the case of $SU(3)$ with $\bar{k}=1$, with $c$ ranging from 0.18 to 0.8. The results for $\mathcal{C}_{1}$ are shown in table 1. Figure 1 displays $\log(\Lambda_{{\rm TGF}}/\Lambda_{\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 sec. 5, but for now we will simply mention that at constant energy scale $\mu=(c\tilde{l})^{-1}$ and fixed $N$, taking $c\rightarrow 0$ is equivalent to taking the large volume limit, in which $\log(\Lambda_{{\rm TGF}}/\Lambda_{\MS})$ should approach the yellow line in the plot.

As for the study of the general dependence on $\hat{\theta}$, we considered a series of coprime values of $\bar{k}$ and (small) $N$ such that $\hat{\theta}$ ranged from 0.14 to 0.5. The full results for $\mathcal{C}_{1}$ are shown in table 2 and fig. 2, in which they are plotted as a function of $\hat{\theta}$ for several values of $c$. We observe that the dependence on $\hat{\theta}$ is rather smooth for the considered values of $\bar{k},N$. A discussion about the $\hat{\theta}$-dependence for larger values of $N$ will be presented in sec. 5.

5 Dependence on the number of colors and the magnetic flux

In this section, we will analyze the dependence of $\lambda(c\tl)$ on the number of colors $N$ and the angular variable $\hat{\theta}=\bar{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. AlvarezGaume:2001tv , in which $N$ is sent to infinity while the torus size is sent to zero in such a way as to keep $\tl$ fixed, and the second is the thermodynamic limit, achieved by simultaneously sending $c$ to zero and $\tl$ to infinity while keeping $c\tl$ fixed. The idea that the infinite volume limit can be attained at $\tl\rightarrow\infty$ by sending either the torus size or the number of colors to infinity is implicit in our construction.

5.1 Singular large $N$ limit and $\hat{\theta}$-dependence

Singular large $N$ limits such as the one described above have been employed in various contexts. In ref. GarciaPerez2015 the non-perturbative running of the $SU(\infty)$ ’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 by $\tl=a\sqrt{N}$, where $a$ denotes the lattice spacing. The continuum limit at fixed $\tl$ 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-comutative gauge theory whose rational adimensional non-commutativity parameter is given precisely by $\hat{\theta}$, a mapping through which the effective torus size $\tl$ corresponds directly to the size of the non-commutative torus in the dual theory. One of the proposals raised in ref. AlvarezGaume:2001tv was to define non-commutative gauge theories at irrational values of $\hat{\theta}$ through a sequence of ordinary $SU(N_{i})$ twisted Yang-Mills theories with increasing number of colors and $\hat{\theta}_{i}=\bar{k}_{i}/N_{i}\rightarrow\hat{\theta}$. In 2+1 dimensions, ref. Perez:2018afi has shown that this is only possible, avoiding tachyonic instabilities, for an uncountable zero-measure set of values of $\hat{\theta}$, such as for instance a sequence of values of $\bar{k}$ and $N$ defined through $\bar{k}_{i}/N_{i}=F_{i-2}/F_{i}$, where $F_{i}$ denotes the $i$th term in the Fibonacci sequence. In that case, instabilities in the large $N$ limit are avoided and the limiting sequence tends to $\hat{\theta}=(3-\sqrt{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}$:

[TABLE]

where the symbol $||x||$ is used to denote the distance from $x$ to the nearest integer Chamizo:2016msz ; Perez:2018afi . 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 Perez:2017jyq .

The limiting procedure to define non-commutative gauge theories at irrational values of the non-commutativity parameter relies on the asumption of continuity in $\hat{\theta}$. The one-loop matching constant $\mathcal{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 ratio $\bar{k}/N$ defining $\hat{\theta}$. While a detailed analysis of the $\hat{\theta}$ dependence is beyond the scope of this paper, we did look at the integrals $I_{1}$ and $I_{2}$ entering the definition of $\mathcal{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 $\mathcal{C}_{1}$ depend on $\hat{\theta}$ 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 $\hat{\theta}$ 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 $\hat{\theta}$ for sequences approaching rational values, in particular for $\bar{k}/N=0$, 1/4, 1/3 and 1/2. A similar dependence on $\hat{\theta}$ 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 Perez:2017jyq , which can be understood in terms of the parameter $Z_{\rm min}$ introduced earlier.

Let us take a look at how the dependence in this $Z_{\rm min}$ quantity enters in the $I_{1}$ contribution to $\mathcal{C}_{1}$. The $\hat{\theta}$-dependent term comes from the function $H(s,u,v,\hat{\theta})$ defined in eq. (101). This contribution is finite in the UV and given by:

[TABLE]

As all terms included in the sum have a non-zero value of $\hat{\theta}\tilde{\epsilon}m$, UV-finiteness is guaranteed. However, in the limit in which this quantity tends to zero, one would retrieve the divergence present in the $\hat{\theta}=0$ term. We will in what follows show that such a limit is approached logarithmically in $Z_{\rm min}$. Let us begin by considering the leading asymptotic behavior for small $x$:

[TABLE]

where $\hat{n}$ denotes the integer closest to $\hat{\theta}\tilde{\epsilon}m$. Integrating over $x$, we get:

[TABLE]

where $Z^{2}(m)=m^{2}\,||\hat{\theta}m||^{2}$. If the argument of the incomplete $\Gamma$ function is small, this goes as:

[TABLE]

The logarithmic dependence in $Z$ is tamed by the exponential damping in ${\hat{c}}m^{2}$, but at small enough ${\hat{c}}$ this suppression disappears, giving rise to the behavior presented in fig. 4(a). This is more clearly seen in fig. 5 where we show the contribution of $I_{1}$ to $\mathcal{C}_{1}$ as a function of $\log Z_{\rm 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_{\rm min}=0.1$. Sequences approaching $\hat{\theta}=0$ in the left plot and $\hat{\theta}=1/2$ in the right one are deep in the region with small $Z_{\rm min}$, where a tiny change in the value of $\hat{\theta}$ translates into a large change in the integral.

As a final remark, we will point out that the value of $Z_{\rm min}$ stays almost constant along the Fibonacci sequence mentioned earlier, meaning that the results of the integrals will 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(\infty)$ running coupling using the reduction techniques employed in ref. GarciaPerez2015 .

5.2 Large volume limit

So far, we have been discussing the dependence of the matching constant $\mathcal{C}_{1}$ on the number of colors and the flux-dependent parameter $\hat{\theta}$ 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 $\mathcal{C}_{1}$ should vanish. As we recall, at fixed value of $t$, $\mathcal{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 $\mathcal{C}_{1}$, and will vanish in the thermodynamic limit provided $\mathcal{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. (59), with $t$ set to $(c\tl)^{2}/8$. As we recall, the dependence on $c$ and $N$ came from:

[TABLE]

where:

[TABLE]

In the infinite volume limit, understood in the sense of $c\rightarrow 0$ at fixed $l_{g}$, one has $F_{0}(0,d)=1$ and therefore:

[TABLE]

leading to a LO term in agreement with the results found in ref. Luscher:2009eq . The leading correction is exponentially suppressed with the square of the volume as:

[TABLE]

If the large $N$ limit (i.e. large $l_{g}$) at fixed $l$ and constant $cl_{g}$ is taken instead, one gets ${\mathcal{A}}(\pi c^{2})=1+{\cal 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 $\mathcal{C}_{1}$ can be written in a compact way as444The regularized expression for $I_{9}$ is slightly different, see app. C.:

[TABLE]

where we used the symbol $\int$ to denote the integrals appearing in eqs. (85)-(96) in a generic manner, including the prefactors multiplying the $\Phi$ function and derivatives when required. The quantity $\hat{H}(s,u,v,\hat{\theta})$ is related to the function $H(s,u,v,\hat{\theta})$ entering the definition of $\Phi$ through:

[TABLE]

and it is given by:

[TABLE]

with:

[TABLE]

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 $\mathcal{C}_{1}$ into $\hat{\theta}-$independent and $\hat{\theta}-$dependent terms, given by:

[TABLE]

where the function $\hat{H}^{\prime}$ is obtained by subtracting the zero modes from $\hat{H}$ after Poisson resummation (see appendix E for the details), and:

[TABLE]

where $z_{\mu}=\epsilon_{\mu\nu}n_{\nu}k/N+ivn_{\mu}/({\hat{c}}N^{2}\alpha u)$, and where $n$ denotes a $d_{t}$-dimensional vector of integers taking values in the intervals $[-N/2,N/2)$ for even values of $N$, and $[-(N-1)/2,(N-1)/2]$ for odd ones. The leading correction to the infinite volume limit is in general driven by the contribution of $I_{TI}^{(0)}$ and $I_{TD}^{(0)}$, and depends on two quantities: ${\hat{c}}\alpha$ and ${\hat{c}}\alpha u/s$.

The simplest case corresponds to integrals $\bar{I}_{1}$, $\bar{I}_{2}$ and $\bar{I}_{4}$, for which both ${\hat{c}}\alpha$ and ${\hat{c}}\alpha u/s$ tend to zero in the ${\hat{c}}\rightarrow 0$ limit in all of the integration range. The leading contribution, derived in appendix E, is given by:

[TABLE]

Integrals for which the infinite volume contribution $I_{i}^{\infty}$ is UV-divergent at $d=4$, such as $\bar{I}_{1}$ and $\bar{I}_{4}$, have a leading correction that goes as $\sim\log(c^{2}N^{2})\exp(-1/(cN)^{2})$. $I_{2}^{\infty}$ 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 fig. 6 the dependence of these integrals on $cN$ for several values of $\bar{k}$ 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 $\mathcal{C}_{1}$ on $cN$ as $cN\rightarrow 0$ is in fact well described by a formula analogous to eq. (168); an example of this for the case of SU(3) is shown in fig. 7, where $\mathcal{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})(\alpha+\beta\log(cN)+\gamma cN+\delta c^{2}N^{2})$. In order to push the calculation of $\mathcal{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 $\bar{I}_{3}$, $\bar{I}_{7}$, $\bar{I}_{10}$, $\bar{I}_{11}$ and $\bar{I}_{12}$, which we were able to compute down to values of $(cN)^{2}\sim 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 $\bar{I}_{1}$, $\bar{I}_{2}$, and $\bar{I}_{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 $\mathcal{C}_{1}$ are considered. We will first analyze the case of $\bar{I}_{3}$ in detail to see how the cancellation takes place, and then generalize it to all other cases. For this integral, in the ${\hat{c}}\rightarrow 0$ limit, ${\hat{c}}\alpha u/s$ goes to zero in the full integration range, and the leading dependence is given by:

[TABLE]

From this expression one can show (see app. E for the details) that the dominant correction in the $cN\rightarrow 0$ limit is:

[TABLE]

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 fig. 8. One can show that, in the infinite volume limit, all remaining integrals $\bar{I}_{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 $\bar{I}_{8}$ and $\bar{I}_{9}$ respectively. Combining eq. (65) with these coefficients, it is easy to show that the total contribution of the leading $(cN)^{2}$ term to $\mathcal{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 fig. 7, we expect other possible power like corrections to cancel out as well when combined to form $\mathcal{C}_{1}$, the final result exponentially decaying towards zero with a leading dependence on $c$ of the form $\sim\log(cN)\exp\left(-(cN)^{-2}\right)/(N^{2}-1)$. A preliminary analysis was performed for the case of $\bar{I}_{3}$, with the quantity $\bar{I}_{3}-I_{0}$ times the factor $(N^{2}-1)\exp((cN)^{-2})$ being shown in fig. 9 as a function of $cN$ for several values of $N$. Each point in that plot was obtained from the exact expression for $\bar{I}_{3}$, and the continuous lines correspond to the approximate expression obtained combining eqs. (164) and  (165). This decomposition is quite useful towards analyzing the $N$ dependence of the integral, and so we displayed each of the two pieces in figs. 10(a) and  10(b) as a function of the appropriate scaling variable $cN$.

The $\hat{\theta}$-independent term is presented in fig. 10(a), multiplied by the factor $(N^{2}-1)\exp((cN)^{-2})$ scaling away most of its $N$ dependence. For $cN\rightarrow 0$, the integral decays exponentially as $\sim\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 $\hat{\theta}$-dependent part is more complicated, as one needs to take into account the dependence on the magnetic flux $k$. The decay of $\bar{I}_{3}$ towards zero is in this case faster than exponential, going as $\sim(cN)^{2}\exp(-(cN)^{-2})$. This is shown in fig. 10(b), where we plotted $\bar{I}_{3TD}^{(0)}$ multiplied by 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 $\hat{\theta}$-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.

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 $\mu=1/(c\tilde{l})$ 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 which $\tl=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 $\Lambda$ 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 $\lambda_{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 keep $\tl$, and hence the renormalized ’t Hooft coupling at scale $\mu=1/(c\tilde{l})$, fixed. This is a singular large $N$ limit in the spirit of those introduced in AlvarezGaume:2001tv , albeit a rather non-standard one since non-planar, $\theta$-dependent diagrams survive the limit as long as $\tl$ 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 length $\tl$ and has a dimensionless non-commutativity parameter given by $\hat{\theta}=\bar{k}/N$. Our analysis also supports the observation, first presented in Chamizo:2016msz ; Perez:2018afi , that the avoidance of tachyonic instablities when taking the singular limit is only possible for a zero-measure, though uncountable, set of values of $\hat{\theta}$. Curiously, one of the successful cases, of limiting parameter $\hat{\theta}=(3-\sqrt{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 Chamizo:2016msz .

The second limit at which we looked was the thermodynamic limit, in which $c$ is sent to zero and $\tilde{l}$ is sent to infinity while keeping the energy scale $\mu$ constant. This leads to the one-loop expression of the ’t Hooft gradient flow coupling at infinite volume Luscher:2010iy . Our results give support to the reduction idea, in the sense that the $SU(\infty)$ coupling in the thermodynamic limit can also be recovered at fixed torus size by sending $N$, and hence $\tl$, to infinity, in which case the limit is approached with $1/N^{2}$ corrections.

Acknowledgments

We would like to thank Antonio González-Arroyo and Alberto Ramos for many valuable discussions on both this topic and related ones. We acknowledge financial support from the MINECO/FEDER grant FPA2015-68541-P and the MINECO Centro de Excelencia Severo Ochoa Programs SEV-2012-0249 and SEV-2016-0597. E.I. Bribián acknowledges support under the FPI grant BES-2015-071791. The numerical computations presented in this paper have been carried out at the IFT Hydra cluster and with computer resources provided by CESGA (Galicia Supercomputing Center).

Appendix 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 Perez:2014sqa and references therein for a review) and in the lattice regularized version of the theory GonzalezArroyo:1982hz ; Luscher:1985wf ; Snippe:1997ru ; Perez:2017jyq . 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. (20), (21), using the irreducible twist given in eq. (11), and the remaining gauge degrees of freedom will be fixed using a generalized covariant gauge of parameter $\xi$ 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:

[TABLE]

where $D_{\mu}\equiv\partial_{\mu}+ig_{0}A_{\mu}$ is the covariant derivative and $c$, $\bar{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. (14), along with an analogous one for the ghost fields:

[TABLE]

where the momenta appearing in these expressions are quantized in units of the effective size $\tl$.

The Feynman rules for the vertices are then obtained from the commutation relations in eq. (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:

[TABLE]

with:

[TABLE]

• •

4-gluon term:

[TABLE]

with:

[TABLE]

• •

Ghost-gluon term:

[TABLE]

These rules can be easily used to derive different quantities, such as the one-loop correction to the propagator. At order $\mathcal{O}(g_{0}^{2})$ and in the Feynman gauge ($\xi=1$), the vacuum polarization tensor can be obtained as shown in ref. Perez:2014sqa :

[TABLE]

Appendix 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 $\mathcal{O}\left(\lambda_{0}^{2}\right)$ contributions to the observable $\left\langle E/N\right\rangle$ in terms of the integrals for this appendix. The contributing ${\mathcal{E}}_{i}$ terms are the following:

[TABLE]

When summing all of the terms contributing to ${\mathcal{E}}^{(1)}(t)$, the $\mathbb{I}_{i}$ terms cancel out, and thus:

[TABLE]

which is the NLO result given in eq. (65). The intervening $I_{i}$ integrals are given in eqs. (85) - (96), in terms of the auxiliary $\Phi(s,u,v,\hat{\theta})$ functions from  (83).

Appendix C Regularization of $I_{9}$

We will present here the full details of the procedure to regularize the integral $I_{9}$, defined in eq. (93), which differs slightly from the general treatment described in sec. 3.2.1. As we recall, the initial integral is split into three terms:

[TABLE]

with the Heaviside function $\theta$ 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}(\Phi^{(0)},t^{\prime})$. The original integral can be rewritten as:

[TABLE]

which, after some algebraic manipulation, becomes:

[TABLE]

where

[TABLE]

This integral ${\cal I}$ is in $d=4-2\epsilon$ dimensions of order $\epsilon$ times zero. The asymptotic behavior at $z=0$ is then obtained by expanding ${{\mathcal{A}}}({\hat{c}}(2t^{\prime}+z))$ around $z=0$, leading to:

[TABLE]

where

[TABLE]

The integral over $z$ presents a pole in $1/\epsilon$, but it is cancelled when multiplied by ${\cal I}$, leading to a final result that is identically zero for $d=4$. This $I_{9}^{\infty}$ 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}(\theta(z-1)\Phi^{(0)},t^{\prime})$ term can be treated in a similar way. One takes the initial expression:

[TABLE]

and rewrites it in a form identical to eq. (187), only with the integral over $z$ restricted to the interval $[1,\infty]$. After some manipulation, the regularized result becomes:

[TABLE]

which is finite in $d=4$ dimensions, and which we were able to evaluate numerically.

Appendix D Numerical implementation of the integration algorithm

In order to perform the numerical computation required to obtain the results presented in sec. 4.2, we prepared a code in C++ to compute the values of the $\Phi(s,u,v,\hat{\theta})$ 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 $\Phi(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:

[TABLE]

which was made finite through the procedure explained in sec. 3.2.1, and where $H$ and $\Phi^{(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:

[TABLE]

where we introduced a generic matrix $X$ to denote either $A_{0}$ or $B$ from eqs. (141) and (144). 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 $\Phi^{\rm 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:

[TABLE]

where we used the shorthands $B_{\theta}\equiv B({\hat{c}}s,{\hat{c}}u,{\hat{c}}v,\hat{\theta})$, $B_{0}\equiv B({\hat{c}}s,{\hat{c}}u,{\hat{c}}v,0)$, and $\tilde{B}\equiv B({\hat{c}}sN^{2},{\hat{c}}u,{\hat{c}}vN,0)$ for clarity. Several of these expressions are redundant:

[TABLE]

allowing us to rewrite the observable $\Phi^{\rm fin}(s,u,v)$ in four equivalent forms:

[TABLE]

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\rightarrow-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 $\Phi^{\rm 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 $\epsilon$ 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 $\Phi^{\rm 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 $\epsilon$ 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\in[0,1]$. The code begins by computing the value of the integrand, which in this case would be the $\Phi$ 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\in[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 $\epsilon$ (the same that we used for the $\Phi$ 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 $\epsilon$ 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 $\Phi$ 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_{fit}=a_{0}-a_{1}e^{-a_{2}(z-z_{max})}$, 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 $\sim 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 $\theta(1-z)$ function being introduced in the integrand (see the end of sec. 3.2.1 and app. 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^{\prime}=1$, but then the Heaviside function became a $\theta(1-y^{\prime}z^{\prime})$ 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^{\prime}$ to be split from $z^{\prime}=1$ onwards into two subintervals $[0,1/z^{\prime}]$ and $[1/z^{\prime},1]$, with the convergence of each side being considered separately.

Due to the procedure we used to determine the convergence of the integrals, for a given integral $\mathcal{I}$, and dubbing the number of integrals to perform $n_{i}$ (single, double or triple), the final error of the integral is:

[TABLE]

This comes from the fact that both the error of the $\Phi$ functions and the convergence criterion for the integrals is given by the same $\epsilon$, so for a single integral:

[TABLE]

Additional integrals simply add extra $1+\epsilon$ 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 $\Delta\mathcal{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 $\Phi$ function returning a value under $10^{-12}$, any inner integral returning any value under $5\times 10^{-12}$ (or $10^{-10}$ in the cases of a few intervals in which using $5\times 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 $\epsilon$ 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 fig. 11 some examples of the profile of the integrand near the aformentioned $(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.

Appendix E The infinite volume and large $N$ limits

In this appendix, we will derive the formulas mentioned in sec. 5, which were used to analyze the $N$ and $\hat{\theta}$ dependence of $\mathcal{C}_{1}$ at NLO in the coupling for the case of a two dimensional twist ($d_{t}=2$). As we recall, the contributions to $\mathcal{C}_{1}$, barring the one from $I_{9}$ which is slightly different, can be written in the form:

[TABLE]

where the notation $\int$ was used to refer generically to the integrals from eqs. (85)-(96), including the prefactors multiplying the $\Phi$ function and derivatives when required. The function $\hat{H}(s,u,v,\hat{\theta})$ was defined through the relation:

[TABLE]

where $\Phi^{\infty}$ was defined in eq. (105) and $H$ is the function given in eq. (101), in terms of which we rewrote the $\Phi$ functions entering the integrals:

[TABLE]

The function $\hat{H}$ can be expressed as:

[TABLE]

with:

[TABLE]

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 $\hat{\theta}$-dependent and $\hat{\theta}$-independent parts. For the latter, Poisson resummation yields:

[TABLE]

In the $\hat{\theta}$-dependent case, on the other hand, we begin by rewriting $m=\hat{m}l_{g}+m^{c}$, with the components of $m^{c}_{\mu}$ 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 to $\hat{m}$ only, leading to:

[TABLE]

where we introduced a $d_{t}$-vector $\chi$ whose components are given by $\chi_{\mu}=||\hat{\theta}\tilde{\epsilon}m^{c}_{\mu}||$, the symbol $||x||$ denoting the distance from $x$ to the nearest integer. Introducing $\chi_{\mu}=n^{c}_{\mu}/l_{g}$ and inverting the relation between $m^{c}$ and $n^{c}$ to write $m^{c}=k\epsilon n^{c}\,($mod $l_{g})$, we obtain:

[TABLE]

The two terms entering $\hat{H}(s,u,v,\hat{\theta})$ and $\hat{H}(s,u,v,0)$ can be rewritten matricially. Recalling the expressions of $A_{0}$ and $B$ from eqs. (141) and (144) we have, in terms of Siegel theta functions, and particularizing to the case of $d_{t}=2$:

[TABLE]

and:

[TABLE]

We will now split the original integral into two pieces, setting $\hat{\theta}=0$ in one part to confine all of the $\hat{\theta}$ dependence to the other one. As we want both of them to be well behaved both in the IR and in the UV, it will be convenient to first isolate the terms corresponding to zero-modes at each step of the calculation, both before and after Poisson resummation. In the original definition of $\hat{H}$, given by eq. (101), the terms with $m=0$ were already subtracted, so we simply need to take away the terms corresponding to $n=0$. Doing so leads to:

[TABLE]

The analogous procedure after Poisson resummation, i.e. subtracting the expressions obtained setting $m=0$ and $n=0$ (separately) in eq. (215), and adding back the $m=n=0$ one, yields:

[TABLE]

where $\hat{H}^{\prime}$ denotes the resulting function after subtracting those zero modes. The same can be done for the term containing ${\mathcal{A}}({\hat{c}}s)$, whose zero mode contribution is given by $-(1-1/N^{2})c^{2}/s^{2}$.

The term $\hat{H}(s,u,v,\hat{\theta})$ containing the $\hat{\theta}$ 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_{\mu}=\tilde{n}_{\mu}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_{\mu}=0$ along periodic directions and $\tilde{n}_{\mu}=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:

[TABLE]

where $z_{\mu}=\epsilon_{\mu\nu}n_{c\nu}k/N+in_{c\mu}v/({\hat{c}}N^{2}\alpha u)$.

We may then rewrite each of the integrals contributing to $\mathcal{C}_{1}$ as the the sum of two components $I=I_{TI}+I_{TD}$, the latter containing all of the $\hat{\theta}$ dependence:

[TABLE]

where:

[TABLE]

and with $n_{c}$ and $z_{\mu}$ as defined above.

From this expression, one can analyze the ${\hat{c}}\rightarrow 0$ limit, whose approach is driven by two variables: ${\hat{c}}\alpha$ and ${\hat{c}}\alpha 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_{TI}^{(0)}$ and $I_{TD}^{(0)}$, and hence the leading order in the ${\hat{c}}\rightarrow 0$ limit for them will be proportional to:

[TABLE]

which approaches zero at least exponentially in the ${\hat{c}}N^{2}\rightarrow 0$ limit, and goes, in the large $N$ limit taken keeping ${\hat{c}}N^{2}$ constant, as $1/N^{2}$. In most cases, the leading contribution in the ${\hat{c}}\rightarrow 0$ limit is hence given by $I_{TI}^{(0)}$ and $I_{TD}^{(0)}$.

The simplest cases are those of $\bar{I}_{1}$, $\bar{I}_{2}$ and $\bar{I}_{4}$, for which both ${\hat{c}}\alpha$ and ${\hat{c}}\alpha u/s$ tend to zero in the whole integral range. Starting from the expressions of $I_{TI}$ and $I_{TD}$, it is easy to derive the leading correction to the large volume limit. In the three cases it is given by:

[TABLE]

All three integrals can be analytically approximated with this, leading to:

[TABLE]

We will now consider the remaining integrals, looking first at the ${\hat{c}}$ dependence of $I_{TI}^{(0)}$ and $I_{TD}^{(0)}$. For $\bar{I}_{3}$, $\bar{I}_{6}$ and $\bar{I}_{10}$, the variable going to zero in the ${\hat{c}}\rightarrow 0$ limit is ${\hat{c}}\alpha u/s$, and the leading dependence is given by:

[TABLE]

whereas for $\bar{I}_{5}$, $\bar{I}_{7}$, $\bar{I}_{11}$, and $\bar{I}_{12}$ the variable going to zero is ${\hat{c}}\alpha$, and we have instead:

[TABLE]

To leading order all these integrals go to zero as $\sim c^{2}$, with a coefficient depending on $N$ that is identical in absolute value for all of them.

We will take a look at $\bar{I}_{3}$ as an illustrative example. The leading contribution in the ${\hat{c}}\rightarrow 0$ limit for this integral is given by:

[TABLE]

which allows us to separate ${\mathcal{A}}$ into two parts, one that depends on $N$ and another that does not:

[TABLE]

Rescaling $z$ to $z^{\prime}={\hat{c}}z$ in the first expression and to $z^{\prime}={\hat{c}}N^{2}z$ in the second, we can decompose the integral into the difference of two pieces $\bar{I}_{3}^{(1)}-\bar{I}_{3}^{(2)}$, which in the $cN\rightarrow 0$ limit become:

[TABLE]

The leading order result in the $cN\rightarrow 0$ limit is thus given by:

[TABLE]

with $a_{1}=-1.76508480122121275$ and for instance $a_{2}(3)=3.59085631503990722$. One can show that all the other integrals are also proportional to $I_{0}$, with the proportionality coefficient being +1 for $i=5,6,7$ and -1 for $i=10,11,12$ respectively. The results for the case of the SU(3) gauge group are displayed on fig. 12, with the red line in the plot showing $I_{0}$ and the remaining continuous lines representing the contribution of $I_{TI}^{(0)}+I_{TD}^{(0)}$. The cases of $\bar{I}_{8}$ and $\bar{I}_{9}$ are shown in the plot as well, which also turn out to be proportional to $I_{0}$ with respective coefficients 4 and -2.