Memory efficient finite volume schemes with twisted boundary conditions

In this paper we explore a finite volume renormalization scheme that combines three main ingredients: a coupling based on the gradient flow, the use of twisted boundary conditions and a particular asymmetric geometry, that for $SU(N)$ gauge theories consists on a hypercubic box of size $l^2 \times (Nl)^2$, a choice motivated by the study of volume independence in large $N$ gauge theories. We argue that this scheme has several advantages that make it particularly suited for precision determinations of the strong coupling, among them translational invariance, an analytic expansion in the coupling and a reduced memory footprint with respect to standard simulations on symmetric lattices, allowing for a more efficient use of current GPU clusters. We test this scheme numerically with a determination of the $\Lambda$ parameter in the $SU(3)$ pure gauge theory. We show that the use of an asymmetric geometry has no significant impact in the size of scaling violations, obtaining a value $\Lambda_{\overline{MS}} \sqrt{8 t_0} =0.603(17)$ in good agreement with the existing literature. The role of topology freezing, that is relevant for the determination of the coupling in this particular scheme and for large $N$ applications, is discussed in detail.


Introduction
Finite volume renormalization schemes [1] are a powerful tool to investigate asymptotically free theories (see [2,3] for a comprehensive review). Asymptotic freedom [4,5] predicts that at high energies the coupling is weak, and perturbation theory becomes a reliable tool to make precise predictions. In Yang Mills (YM) theories or QCD, the approach to this perturbative regime is logarithmic. Due to this logarithmic running of the coupling, the strongly coupled regime and the perturbative regime are separated by large energy scales. Accommodating these disparate scales on a single lattice simulation is challenging, and requires compromises (see [6] for a review on the consequences in the context of the extraction of the strong coupling). Finite volume renormalization schemes solve this multiscale problem by linking the physical volume of the system with the renormalization scale. A single lattice simulation only allows to study a limited range of scales, but by using a recursive procedure called finite size scaling, simulations with different physical sizes allow to cover large energy scales. Using this technique, the running coupling has been determined in QCD with N f = 2, 3 [7][8][9] and 4 [10] flavors as well as in the pure gauge theory [11][12][13].
In this work we will explore numerically a particular finite volume renormalization scheme, first introduced in [14][15][16] based on three ingredients. First, we choose a nonperturbative coupling based on the gradient flow [17][18][19][20] (see also [21] for a review). Second, we choose twisted boundary conditions [22][23][24] for the gauge fields. Finally, we make use of a particular geometry, namely, for the case of SU (N ) gauge theories we choose a hyper-cubic box of dimensions l 2 ×l 2 withl = l N .
The logic behind these choices is rooted in the study of the large N limit of YM theories and in the ideas of volume reduction with twisted boundary conditions [25][26][27] (see [28,29] for recent reviews). With our choice of geometry (l 2 ×(Nl) 2 ) and our particular choice of twisted boundary conditions, perturbation theory indicates that the effective volume of the torus is (Nl) 4 . This equivalence between group degrees of freedom and spatial degrees of freedom (that becomes exact in the large N limit), lies behind the ideas of volume reduction and non-commutative geometry.
In this work we will will argue that this particular choice of scheme has several advantages. The running coupling in the pure gauge theory has recently seen a revived interest in relation of the quest for precision in the LHC. It has been shown that precise values of the parameter of the pure gauge theory can be used, via a non-perturbative matching between QCD and the pure gauge theory using heavy quarks [30], into a precise value for the strong coupling, a key quantity for phenomenology in high energy physics. Due to the fact that the pure gauge theory is numerically more tractable than QCD this approach is probably the best option to substantially reduce the uncertainty in the strong coupling. Still, very precise determinations of the running coupling in the pure gauge theory are needed, that poses its own challenges. Compared with the more customary finite volume renormalization schemes based on Schrödinger functional (SF) [31,32] or mixed open-SF boundary conditions [33], twisted boundary conditions [34,35] preserve the invariance under translations and are therefore free of the linear O(a) cutoff effects present in other schemes. The use of periodic boundary conditions [36,37] lead to coupling definitions that are non-analytic, making cumbersome the extraction of the parameter. This is not the case for twisted boundary conditions, which are particularly suited for perturbative calculations. Moreover, the relation between the parameter in this scheme and in the MS scheme is known analytically [14]. Finally, due to the use of an asymmetrical torus, the memory footprint required to simulate this scheme in a computer is reduced. The comput-ing power stays roughly the same because with this geometry the observables show a larger variance [38], but this scheme allows for a better usage of current GPU clusters.
In this work we will determine the parameter in the pure gauge theory. We will check that the use of this particular geometry has no significant impact in the size of scaling violations, and will explore in detail the role of topology freezing [39,40]. In summary we expect that this particular scheme will play an important role in future precise determinations of the strong coupling, and in the investigation of the large N limit of gauge theories.
Let us now comment a bit on our choice of coupling and geometry. Definitions of the gauge coupling based on the gradient flow are nowadays customary since they can be computed on the lattice to high precision. The gradient flow is a renormalization procedure based on replacing the original gauge fields A μ (x) by a new set of smooth, flow time dependent, fields B μ (x, t), in terms of which gauge invariant composite observables are automatically renormalized quantities for t > 0 [41]. These new fields are driven by the so-called flow time equations: where D μ and G μν stand for the covariant derivative and field strength tensor of the flow fields. The effect of the flow is to smooth the gauge fields on a sphere of radius √ 8t. Renormalized couplings based on the gradient flow can thus be trivially introduced, as one only needs to find dimensionless observables that depend on a scale. Starting from E(t), the energy density of the flow fields: we can use the flow time t, with dimensions of length squared, to construct the dimensionless combination t 2 E(t) . Nevertheless, this dimensionless combination depends on a scale μ = 1/ √ 8t, which allows to define a renormalized 't Hooft coupling (λ = g 2 N ) through the relation: where N stands for a normalization constant introduced to ensure that λ(μ) = λ MS (μ) + · · · at leading order in perturbation theory. In the context of finite-size scaling, it is customary to relate the renormalization scale to the size of the finite volume in which the gauge theory lives. In the standard set-up, on a torus of size l 4 , one fixes the flow time via the relation √ 8t = cl, with c standing for an arbitrary, smaller that 1, constant defining the scheme.
Concerning the geometry, the twisted gradient flow (TGF) scheme is defined by introducing SU (N ) Yang-Mills theories on an asymmetrical hypercubic box of size l 2 × (Nl) 2 .
The plane with two short directions, taken to be the (1,2) plane, is endowed with twisted boundary conditions, i.e.
In the remaining directions the gauge field is periodic with periodl = Nl. The twist matrices ν , for ν = 1, 2, satisfy the consistency condition: with k and N coprime integers. For the concrete case of SU (3) studied in this paper, the choice of twist is unique and given by k = 1. However, for arbitrary N the value of the twist k has to be scaled with the range of the gauge group and is best selected so as to avoid the appearance of tachyonic instabilities in the large N limit [27,42]. As already mentioned, perturbation theory indicates that finite volume effects for this choice of twist are controlled bỹ l [14,43,44]. Under this assumption, it becomes natural to use the effective sizel, rather that l, to set the scale of the running coupling [45]. In addition, and following the prescription introduced in [40], the coupling will be defined within the sector of configurations with zero topological charge. This choice aims at circumventing the problem of topology freezing on the lattice that will be discussed in some detail in Sect. 4. With this, the coupling is defined as: where δ Q stands for a δ-function that restricts the path integral to configurations with zero topological charge Q, and where θ 3 (0, i x) stands for the Jacobi θ 3 function: Throughout this work we make use of the coupling scheme corresponding to c = 0.3. The paper is organized as follows. Section 2 summarizes the strategy to determine the parameter in terms of some typical hadronic scale. Section 3 compiles our numerical results. In Sect. 3.1 we introduce our numerical set-up and the lattice definition of the TGF coupling. Section 3.2 discusses the determination of the step scaling function and the continuum extrapolation. In Sect. 3.3 we determine the parameter in the TGF scheme in terms of a hadronic scale μ had , and discuss the matching with the SF and MS schemes.
Our final result for MS √ 8t 0 is presented in Sect. 3.4. Finally, Sect. 4 discusses the role of topology freezing. We conclude in Sect. 5 with a summary of results. A few appendices are included providing raw data for our lattice results, details on the matching to the SF scheme and a discussion of the dilute gas approximation used in Sect. 4.

The Lambda parameter
The starting point for the determination of the parameter is the RG equation in a certain scheme (labeled s) .
The β-function has a perturbative expansion with the first two coefficients (4π) 2 b 0 = 11/3 and (4π) 4 b 1 = 34/3 being universal (i.e. scheme independent). The first order differential equation (2.1) can be integrated exactly where s is an integration constant (the -parameter). As such, it is a renormalization group invariant (i.e. scale independent), but depends on the choice of renormalization scheme. This scheme dependence can be computed exactly via a one-loop computation. In particular, the quotient of parameters in the MS scheme and the TGF one (see Sect. 1 for a precise definition of the latter) is given by [14] log TGF MS = 3 22 where C 1 (c = 0.3) = 0.508(4) for gauge group SU (3) and twist k = 1.
In order to determine the parameter, we use the exact relation Eq.
, (2.5) where μ had is a typical hadronic scale and μ pt μ had is a high energy scale where perturbation theory is applicable. Lattice field theory allows to determine non-perturbatively the running of the coupling in some suitable schemes (like the TGF used here), and therefore the second factor of the right hand side of Eq. (2.5). In particular the technique of finite size scaling [1] allows very precise determinations. We postpone a detailed discussion for the next section, and focus here on the evaluation of the first factor, that requires the use of perturbation theory. If we define (2.7) The corrections decrease as powers of λ s (μ pt ), and therefore logarithmically (i.e. slowly) with the scale μ pt . The extraction of the s parameter has to be performed by taking the limit .
the step scaling function allows to determine precisely the last factor in Eq. (2.8). Once the function σ (u) is known, one defines u 0 = λ s (μ had ), and recursively determines Each application of the step scaling function changes the renormalization scale by a factor 2, so that u k = λ s (2 k μ had ) and After a reasonable number of steps (O(10)), large energy scales have been achieved, where perturbative corrections are expected to be small and one can explore the limit of Eq. (2.8) by taking λ s (μ pt ) = u k and μ pt = 2 k μ had μ had . Crucially, large energy scales can be covered using this procedure. In the next section we discuss the numerical determination of the step scaling function in the TGF set-up.

Numerical set-up
We simulate pure SU (3) Yang-Mills theory on asymmetric lattices of size L 2 ×L 2 , withL = N L ≡ 3L, corresponding to a torus in the continuum of size l 2 ×l 2 , withl = aL and a the lattice spacing. Twisted boundary conditions are implemented by taking as lattice action with plaquettes P μν (n) given in terms of the SU (3) link variables U μ (n) as: and where b = 1/λ 0 stands for the inverse of the lattice bare 't Hooft coupling, and Z μν (n) is set to one for all plaquettes except for the ones with coordinates x 1 = x 2 = 1, where: We use a hybrid over-relaxation algorithm [46], that combines a single heat-bath sweep (HB) [47,48], followed bỹ L over-relaxation sweeps (OR) [49]. Since the measurement of flow quantities are computationally expensive, we per-formL hybrid sweeps between measurements, which produces uncorrelated measurements basically for all our coupling values (nevertheless, the small autocorrelations are still taken into account). The simulated values of the bare coupling and lattice sizes as well as the total number of configurations attained in each case is reported in Table 6 in Appendix A. We have used lattices withL = 12, 18, 24, 36 and 48, in a range of values of b ∈ [0.33, 1.12]. Data analysis has been performed using two different approaches. One based on resampling techniques and the other using the method [50] together with automatic differentiation [51]. Both approaches give similar results in all cases.
In this study, we have used the so-called Wilson flow discretization of the continuum flow equation, and evaluated the energy density by means of the clover discretization given by: with P μν (n, t) the plaquette evaluated in terms of the flowed links at flow time t.
It is well known that simulations at very fine lattice spacings loose ergodicity [39]. In many of our simulations this is not problematic, since they are performed in small physical volumes, where the contributions of the non-trivial topological sectors are highly suppressed. But at physical volumes l ∼ 1 fm, critical slowing down can severely affect the results of our studies [40]. For this reason, we define the coupling in the zero topological sector using: where the topological charge Q on the lattice is determined through the expression: evaluated at flow time √ 8t = cl. Since the topological charge given by this expression is not an integer, we have set: It has to be emphasized that the renormalization of the coupling with and without projecting to the trivial topological sector is different. Note that the projection to the trivial sector is performed with the flowed definition of the topological charge (c.f. Eq. (3.7)), that has a well defined continuum limit. Therefore no additional divergences will arise with this definition. Moreover since both coupling definitions have the same perturbative expansion to all orders in the continuum, we conclude that the difference has to be non-perturbative finite terms: just a different choice of scheme.
In addition, and in order to eliminate the leading order lattice artefacts in perturbation theory, the lattice normalization factor has been used: (3.9) whereq μ = 2 sin(q μ /2) stands for the lattice momentum, with q μ = 2π n μ /L, with n μ = 0, . . . ,L − 1, and with the prime in the sum denoting the exclusion of momenta with both components in the twisted plane satisfyingLq i ∝ 2N π = 6π . The values for the TGF coupling measured for c = 0.3 are collected in Table 7 in Appendix B. At a given value of b, the coupling is measured on two lattices with extents given bỹ L and 2L, which allows to determine the lattice step scaling function (u,L) forL = 12, 18, 24.

Determination of the continuum step scaling function
Once the lattice step scaling function (u,L) is determined, one needs to extrapolate it to the continuum. In our set-up, leading cutoff effects are O(a 2 ). Even though this leading behavior receives logarithmic corrections [52], our range of lattice sizes span a factor two in lattice spacings, so we expect these logarithmic corrections to be small. Since the aim of this work is not to provide a very precise result, but else to explore the viability of the TGF scheme, we will rely on simple linear extrapolations. If data at fixed values of the coupling u i (i = 1, . . . , N points ) and their corresponding values of the lattice step scaling function i (L) ≡ (u i ,L) exist for different lattice sizesL, we can extrapolate them as: Alternatively, we have tried to extrapolate The difference between the two approaches is clearly an 1/L 4 term, and in our checks it turns out to be negligible. These extrapolations at fixed values of u i result in the corresponding values σ i = σ (u i ) of the step scaling function. This approach is described in detail in Sect. 3.2.1. It is also convenient to find out a suitable parametrization for the continuum function σ (u), in order to obtain the sequence of couplings Eq. (2.12). Reasonable parametrizations for the functional form σ (u) can be found by using perturbation theory as a guide. In particular, we will use the simple parametrization to fit our data. Note that although this functional form imposes the correct perturbative behavior (for u → 0), after determining the fit parameters p k using our non perturbative data, we expect it provides a description of the nonperturbative step scaling function σ (u) in our range of couplings. Alternatively, being σ (u) a one-to-one function in our domain of interest, one can use the inverse relation to parameterize the inverse step scaling function σ −1 entering in Eq. (2.12), i.e.
The two parametrizations should give equivalent results. Finally we can combine both, the continuum limit and the parametrization of the step scaling function, in a single global fit. This approach does not require the coupling to be perfectly tuned to constant values of u, and amounts to fit the data for the lattice step scaling function to the functional and similarly for the parametrization of the lattice inverse function Details of this global approach are presented in Sect. 3.2.2.

Continuum limit on a u-by-u basis
We have tuned the bare coupling b on theL = 12, 18 and 24 lattices to attain a few targeted values of the renormalized coupling u (see Table 8 in Appendix C). Once this is achieved, the renormalized coupling is measured on the double size lattices at the same value of the bare coupling. Since the tuning is not perfect, there is a small mismatch in the tuning of the target couplings that can be easily corrected for by slightly shifting the resulting couplings to a constant value of u. By using the one-loop perturbative relation the values of the lattice step scaling function (u,L) can be shifted to a target value of the coupling u tg using the relation We note that our data set is already very well tuned, and therefore the additional shift introduced by this procedure is well below the statistical accuracy. The same one-loop relation can be used to propagate the error of u into an error of the step scaling function. Being precise, we add in quadratures to the statistical error of (u,L). This procedure leads to a set of data at constant value of the coupling u for different lattice spacings. The result can be seen in Table 1, while the raw data is available in Table 8 in Appendix C.
The continuum limit of the step scaling function is obtained by performing a linear extrapolation in 1/L 2 , as illustrated in Fig. 1. The continuum values of σ (u) resulting from these fits are given in Table 1. The advantage of this procedure vs the global fit presented below is that no functional dependence in u is assumed for the terms representing the lattice artefacts (cf. Eqs. (3.11) and (3.14)).
The resulting pairs (u, σ ) are then fitted to the series expansions given by Eqs. (3.20)

Continuum limit from a global fit
As an alternative to the u-by-u fit, the continuum limit can be obtained from a combined parametrization of the continuum step scaling function and a continuum extrapolation -see e.g. Eq. (3.14). Table 1 Values of the lattice step scaling function (u tg ,L) at fixed targeted values of the coupling u tg and for different values ofL. We also report on the continuum extrapolation σ (u tg )  (10) 13.91650 32.01 (21) 33.21 (27) 34.00 (25) 34.52 (29) Fig. 1 We display the continuum extrapolation of (u)/u, performed at a few selected values of u From the last section, we know that the continuum step scaling function is well parametrized using three parameters in addition to the two universal terms (i.e. N c = 4). Cutoff effects are well described as long as N l ≥ 3, with little differences in the continuum results even if the number of parameters is increased up to N l = 8. The raw data used in these fits is available in Appendix B. We have fitted the data for 1/U(σ,L) vs σ in the range σ ∈ [1.5, 17] with the particular choice N c = N l = 4; the resulting continuum fit parameters are given by: The result of this fitting procedure is displayed in Fig. 2a where the raw data for 1/U(σ,L) is represented vs σ and compared with the continuum extrapolation resulting from this global fit which is given by the pink band displayed in the plot. We finally present a comparison between the various continuum extrapolations in Fig. 2b. The orange points come from the direct extrapolation at fixed targeted values of u given in Table 1. The blue and pink bands correspond respectively to the u-by-u and global fits described previously, showing a good agreement in all the range of fitted values. For comparison, the perturbative predictions at one and two-loop order are also shown in the figure.

Determination of TGF /μ had
The determination of TGF /μ had relies on the use of Eq. (2.8), evaluated in terms of a sequence of couplings connecting λ(μ had ) and λ(μ pt ). We have performed this sequence in two steps: an exact step connecting the coupling at scales μ had and 2μ had , followed by a step-scaling sequence, starting at λ(2μ had ), to be determined by the con-tinuum step scaling functions obtained in the previous sections. Let us now describe these steps in detail.

Determination of λ(μ had )
We fix the coupling λ(2μ had ) ≡ 13.9164955 at values of the lattice spacing that correspond toL = 12, 18, 24. The corresponding values of the lattice step scaling function are given in Table 1, and extrapolate to a continuum value λ(μ had ) = 34.52 (29). Determination of MS /(2μ had ) Using the fits of the continuum step scaling function, we determine a sequence of couplings λ n = λ(2 n+1 μ had ). The differences between the global and the u-by-u fits in this step are well below our statistical uncertainties. On the other hand, there are several possibilities to determine the ratio MS /(2μ had ): TGF: One can use the universal coefficients of the βfunction and the known ratio MS / TGF to determine MS /(2μ had ) (cf. Eq. (2.8)). MS: One can use the value of the ratio MS / TGF to determine the coupling in the MS scheme as follows (3.23) Then, using the known 5-loop β-function [53] in the MS scheme, one can determine the ratio MS /(2μ had ). SF: Very precise data for the SF coupling is available in reference [13]. These computations were performed at the same values of the bare coupling b as our runs but on lattices of sizeL/3. This allows to convert non-perturbatively the values of λ(2 n+1 μ had ) to values of λ SF (0.9 × 2 n+1 μ had ). Details of this matching are available in Appendix D, here it is enough to say that this matching is only possible in the region 9 λ 3. After this non-perturbative conversion, the known 3loop β-function in the SF scheme [54][55][56] can be used to determine the ratio MS /(2μ had ).
This procedure is similar to the approach taken in [38], and we encourage the reader interested in the details to consult that reference.
The results for the ratio MS /(2μ had ) using the different procedures described above can be seen in Table 2 and in Fig. 3. Methods labeled TGF and MS carry a perturbative uncertainty O(λ(2 n+1 μ had )) (right panel of Fig. 3), while the SF method carries an uncertainty O(λ 2 (2 n+1 μ had )) (left panel). Figure 3 shows that these perturbative corrections are significant in the methods labeled TGF and MS. Nevertheless the picture shows a nice agreement between all analysis methods once the limit λ → 0 is taken.
Finally, we quote as central value of the analysis the one labeled as SF, since it has the smallest perturbative corrections and is compatible within errors with the other two determinations. We report as final numerical value: In order to determine the dimensionless combination MS √ 8t 0 , we have to put together our previous results for the ratio MS /μ had with a determination of μ had × √ 8t 0 . This quantity is computed as the continuum extrapolation of where aμ had is determined by the lattice size and bare coupling for which the renormalized coupling takes the value λ(μ had ) = 34.52(29) (see Table 1). Values for √ 8t 0 /a for our choice of Wilson gauge action are available in the literature (see Table 3).
Two steps were therefore required. First, for each value of L we determined the values of the bare coupling leading to a renormalized coupling λ(μ had ). For this purpose, we have fitted the dependence of the coupling on b, at fixedL, with a Padé fit of the form: (3.26)  Table 3 Results for t 0 /a 2 and r 0 /a for different values of b. For the case of t 0 /a 2 the relevant references are † [20], ‡ [57], [58]. For r 0 /a, instead, the results are from † † [59] and [58]. The data labelled § is obtained from the values of r c /a of [60] together with r c /r 0 = 0.5133(24) which gives the quoted values of r 0 /a Fits are done to 1/λ − b, as illustrated in Fig. 4. After this, one has to determine a/ √ 8t 0 at those values of the bare coupling. Following Ref. [13], this is done by using the results quoted in Refs. [20,57,58]. All the values for this quantity given in table 4 of Ref. [13] are used to fit log(t 0 /a 2 ) to a polynomial in b of degree 5 in the range b ∈ [0.33, 0.39]. This fit is then used to determine a/ √ 8t 0 at the specific values of the bare coupling corresponding to λ(μ had ) on each lattice. The results obtained from this procedure on lattices withL = 12, 18, 24, 36 and 48 are given in Table 4. Figure 5a shows the continuum extrapolation of μ had √ 8t 0 obtained using this strategy. We include datapoints corresponding to theL = 18, 24, 36 and 48 lattices. The final continuum extrapolated value is given by: This number is in agreement within errors with a more precise determination obtained in a recent work by one of the present authors [13]: √ 8t 0 × M S = 0.6227(98). The strategy to obtain the parameter in that work is similar to the one presented here but is based on the use of a different gradient flow scheme with SF instead of twisted boundary conditions.
A more detailed comparison with the existing literature is better done in terms of the scale r 0 [67], which, although less precise that t 0 , is used in most of the determinations of the parameter up to date. An analysis completely analogous to the one leading to Eq. (3.28), allows to determine r 0 × MS . The values of r 0 /a employed in the extrapolation are shown in Table 3. The values of the lattice spacing a/r 0 corresponding to the renormalization condition λ(μ had ) are shown in Table 4 and the continuum extrapolation is displayed in Fig. 5b, leading to a value of μ had ×r 0 = 1.453 (19).
Our final result for the parameter in units of r 0 is: (20).
(3.29) Figure 6 shows a comparison of this result with the values included in the FLAG average [61]. With the present accuracy, our result comes out compatible with both the FLAG average r 0 × MS = 0.615 (18), and the result by Dalla Brida and Ramos: r 0 × MS = 0.660(11) [13]. We believe that our particular scheme can be efficiently used to pin down the error in the determination of the parameter and clarify the tension among the existing values in the literature.

The effect of topology
As explained in the introduction, our definition of the gradient flow coupling follows the prescription adopted in [40] which amounts to evaluate the coupling in the sector of trivial topology. This choice can be seen as a particular renormalization scheme, introduced to circumvent the well known problem of topology freezing and critical slowing down [68]. A neat example of this effect is given in Fig. 7 where we display the value of the topological charge as a function of Monte Carlo time for lattices withL = 24, 36, 48 and a physical size corresponding to λ ≈ λ(μ had ) (l ∼ 1.1 fm). It is clear that our statistics in the lattices with finer lattice spacing is not enough to correctly sample the topological charge. Furthermore, as Table 5 illustrates, the correlation between the coupling and the topological charge in these ensembles is very strong: values of the coupling averaged over the Q = 0 or |Q| = 1 sectors differ by as much as 30%. Clearly, when topology is poorly sampled, any attempt to compute the coupling in this volume regime by averaging over all topological sectors would lead to biased results.
In the remaining of this section we will argue that even when topology is frozen, topological charge fluctuations are correctly sampled in the trivial topology sector (see [69] for a study in the Schwinger model). We will also make an attempt to explain the observed dependence of the coupling on topology in semiclassical terms.   [61]. These are: ALPHA 98 [62], QCDSF/UKQCD 05 [63], Brambilla [64], Flow 15 [65], Kitazawa 16 [66], Ishikawa 17 [12]. We also include the recent work Dalla Brida, Ramos 19 [13] in the comparison

Topological charge fluctuations in the sector of trivial topology
In order to determine whether local fluctuations of topological charge are correctly sampled in the trivial topology sector, we have computed the two-point function of the operator: obtained by integrating the topological charge density, smeared over a flow time t, over three directions including the two short ones. On the lattice, this operator has been determined using the field theoretical definition of the charge density, cf. Eq. (3.5), Tr G cl μν (n, n 0 , t) G cl μν (n, n 0 , t) , (4.2) evaluated at flow time √ 8t = 0.3l. In terms of this operator, we computed the correlator: defined over the sector with trivial topology. Figure 8a shows the time dependence of the correlator evaluated on the Monte Carlo ensembles displayed in Fig. 7. Even though topology is frozen on the largest lattice with L = 48, we observe no significant deviation from scaling; the structure of the charge correlation is the same on the three lattices. These good scaling properties are independent of the physical size of the box and preserved at smaller volumes, an example forl ∼ 0.66 fm is displayed in Fig. 8b. These results indicate that local fluctuations of topological charge are correctly implemented in the trivial topology sector, even when the total topological charge is frozen to zero in the Monte Carlo ensemble.

Semiclassical picture
The step scaling sequence interpolates between the perturbative domain, where non-trivial topological charge is suppressed, and the non-perturbative, strongly coupled regime. In between, in the small to intermediate volume domain, one expects semiclassical arguments based on instantons to hold. In this subsection those arguments will be used in an attempt to reproduce some of the features observed in the simulations, in particular the correlation between coupling and topological charge. By no means should this discussion be interpreted as a rigorous proof but just as a compelling qualitative picture of the role that instantons may play in this context.
In Appendix E we discuss how to apply the dilute instanton gas picture to compute the leading semiclassical contribution to the TGF coupling in a fixed topological charge sector. The derivation is based on the observation that the flow removes short distance fluctuations but preserves classical solutions of the Euclidean equations of motion. This includes instantons although not instanton-anti-instanton pairs, the latter will only be preserved at typical distances larger than √ 8t = cl. The leading contribution to the TGF coupling in the semiclassical approximation can be estimated taking into account that: This expression can be combined with the fact that the classical Euclidean action of a configuration with well separated instantons and anti-instantons can be approximated by S 0 (n +n), with n(n) the number of instantons(antiinstantons) and S 0 the one-instanton action. Therefore, at lowest order in the dilute gas approximation, the energy density receives a (semi-)classical contribution that is proportional to the average number of instantons plus antiinstantons, i.e., This contribution adds to the perturbative one. One can as well evaluate in this approximation the average value of any observable, in particular n +n , restricted to the sector of topological chargeQ by using: with Q the total topological charge, in analogy to what we did for the definition of the coupling in the zero topological sector in Eq. (1.7). In order to determine the average number of instantons in the dilute gas approximation, we have considered two possibilities: a dilute gas of ordinary instantons of action S 0 = 8π 2 , and one made out of fractional instantons [70] with instanton action S 0 = 8π 2 /N . The derivation follows the standard treatment, see i.e. [71]. In this context, the semiclassical picture based on the less standard case of fractional instantons follows the ideas put forward by González-Arroyo and collaborators in a series of works [72][73][74][75][76], see also [77][78][79].
Details of the calculation are provided in Appendix E, here we will only summarize the main final formulas. The contribution to the coupling derived in this approximation, when restricted to the sector of topological charge Q, reads: with ν = Q for ordinary instantons, and ν = N Q for fractional ones, and where A(x) = A(x)/x 2 , coming from the coupling normalization, cf. Eq. (1.7). In this expression, I ν (x) stands for the modified Bessel function of the first kind, and x = 2RV , with R the probability to generate one instanton, or anti-instanton, per unit volume. A way to determine R from the numerical simulations is to measure the topological susceptibility. In the dilute gas approximation and for the case of ordinary instantons, the relation between R and χ is simply given by: χ = 2R, while for fractional SU (3) instantons one obtains instead: The topological susceptibility is defined formally in the continuum through the correlator of the topological charge density q(x): A proper definition requires to regulate the short-distance singularity in this expression, and this can be attained by evaluating the charge density at non-zero flow time. On the lattice, in the absence of topology freezing, the topological susceptibility can be computed by measuring: with Q the total topological charge, determined from the field theoretical definition Eq. (3.7), evaluated on the flowed gauge fields. An illustration of the volume dependence of the topological susceptibility is provided in Fig. 9 where we display t 2 0 χ as a function ofl/ √ 8t 0 . The results shown correspond to simulations with lattice spacing a > 0.03 fm (we refrain to give errors in this plot since our determination of the susceptibility may be affected by freezing, in particular for values of a < 0.05 fm). The onset of the susceptibility is very abrupt and takes place for values ofl/ √ 8t 0 ∈ [1, 3], overlapping with the region covered by the step scaling sequence. After that, χ stabilizes to the constant large volume value, given in SU (3) by t 2 0 χ = 6.67(7) × 10 −4 [80], corresponding in the figure to the orange horizontal band.
By fixing x in terms of the measured susceptibility, one can evaluate the dilute gas contribution to δλ(c) = λ Q=1 (c) − λ Q=0 (c) and compare the result to the values of δλ(c) determined in the simulations. This has been done for the set of Fig. 9 Volume dependence of the topological susceptibility evaluated from flowed fields at 8t = (0.3l) 2 . The orange horizontal band corresponds to the value of the, continuum extrapolated, large-volume SU (3) susceptibility t 2 0 χ = 6.67(7) × 10 −4 [80]. The red line gives the dilute gas estimate discussed in Sect. 4.2, cf. Eq. (4.11) with C = 880 simulations with β ∈ [6, 6.8] where topology seems to be properly sampled and the number of Q = 1 configurations is sufficiently large. In addition, the calculation has been done for various values of the flow time: 8t = (cl) 2 , with c = 0.2, 0.3 and 0.4. We observe that the quantity A(π c 2 )δλ(c) is quite independent of the value of c. Our results are displayed as a function of χ V in Fig. 10 and compared with the dilute gas predictions for ordinary and fractional instantons, derived from Eq. (4.7). Surprisingly, the dependence predicted by the fractional instanton formula matches remarkably well our results in a large range of volumes. We stress that the lines displayed in the figure do not have any free parameter, they are directly given by Eq. (4.7) with the input parameter x fixed in terms of the measured value of the susceptibility.
In the dilute gas approximation one can as well obtain a prediction for the parameter x = 2RV in terms of the renormalized coupling. As explained in Appendix E, the result for fractional instantons at one loop order reads: where ρ stands for the characteristic fractional instanton size, proportional tol in the small volume regime, and C for a normalization constant, unfortunately unknown. Disregarding this fact, we have used this functional form to get a qualitative idea of the volume dependence of the topological susceptibility in this approximation. The calculation requires to input the renormalized coupling constant at a scale ρ ∝l; differences in the scheme dependence translate at this order in a change of value of C that depends on the ratio of parameters. We have tested the formula setting ρ = cl and using the coupling in the TGF scheme with c = 0.2, 0.3 and To summarize, although the results presented in this section have to be taken with a grain of salt, as they are not based on a rigurous analytic calculation, we believe they provide compelling indication that some of the observed effects, in particular the correlation between topological charge and coupling, can be understood in semiclassical terms.

Conclusions
In this work we have investigated the running coupling using finite size scaling techniques in a scheme with twisted boundary conditions. We have used a particular choice of geometry, motivated by the study of the large N limit of SU (N ) Yang Mills theories. In this scheme the physical size of the four dimensional volume is asymmetric, with two directions being a factor N larger than the other two.
We have focused our investigations in the group SU (3). For this case there are many results available in the literature, which allows for a more stringent comparison. Our determination of the MS parameter in units of the flow scale √ 8t 0 or the Sommer scale r 0 shows a good agreement with the literature. The precision of our computation does not match the one of the most precise determinations, but the systematic associated with the continuum extrapolation is taken into account seriously and our simulations involve large lattices. In addition, the matching with perturbation theory is performed in several ways, using a non-perturbative running that covers large energy scales.
In this respect, it is important to note that the determination of MS is accomplished in two different ways. First, we perform a direct determination in our finite volume scheme with twisted boundary conditions. In this case the ratio MS / TGF is a key piece in the computation. Second, we perform a non-perturbative matching with the SF coupling, using available data from the literature. The two determinations show a good agreement, but they exhibit a subtle point. The extraction of the parameter at a certain energy scale μ pt has perturbative corrections O(λ k (μ pt )) (with k = 1 in the case of the direct extraction, and k = 2 for the case of the non-perturbative matching with the SF). A crucial point to obtain compatible results is to remove these subleading perturbative corrections by performing an extrapolation λ(μ pt ) → 0. At our level of precision, and specially for the case of the direct extractions, perturbative corrections are significant, even at energy scales ridiculously large.
All in all, this work shows two promising directions for future research. First, extracting the N dependence of the -parameter in this particular scheme is computationally cheaper than a brute force approach. Second, the particular choice of geometry shows a smaller memory footprint than usual simulations with symmetrical lat-tices. This smaller memory footprint implies an increase in the ratio computing over memory transfer, and therefore can be executed more efficiently in current hardware accelerators. Here it is important to stress that this pure gauge case can be non-perturbatively matched to QCD using heavy quarks [30] and has therefore direct phenomenological relevance for the extraction of the strong coupling in QCD.
The authors hope to investigate these points in detail in future works.

Appendix C: Continuum extrapolation for the step scaling function
See Table 8.  In order to match the values of the coupling λ TGF (μ) to the SF scheme, we use a set of simulations performed at the same values of the bare coupling but on a symmetric lattice with size L SF =L/3 with SF boundary conditions and an abelian background field [31]. The values of the SF coupling in these simulations can be checked in reference [13]. 1 We perform a fit, discarding the SF data withL = 18, of the form where c k are the fit coefficients and The result of the fit can be seen in Fig. 11. This fit can only be trusted for 3 < λ TGF < 6, since data on the SF scheme at higher energies is not available in the literature. 1 Reference [13] gives the values of the SU (3) Yang-Mills couplingḡ 2 SF related to the 't Hooft coupling used in our work through the standard relation: λ SF = 3ḡ 2 SF .

Appendix E: Dilute gas approximation
In this appendix we present the derivation of the dependence of the coupling on the topological charge in the dilute gas approximation. The starting point is the continuum expression for the TGF coupling: with V the volume of the asymmetric box, one gets: We will use that: In a configuration with well separated instantons and antiinstantons, the classical Euclidean action can be approximated by S 0 (n + n), with S 0 the one-instanton action and n(n) the number of instantons(anti-instantons). Therefore, one can compute the contribution to the expectation value of the energy density at leading order in this approximation

Fig. 11
Non perturbative matching between the SF and the TGF couplings. We use the available data for the SF coupling in [13]. In the fit we discard the data withL = 18, but these are plotted to show the consistency of the fit by estimating the average number of instantons plus antiinstantons in the flowed ensemble: Here we rely on the assumption that the flow removes short distance fluctuations but preserves instantons, note though that this is not the case for instanton-anti-instanton pairs which will only be preserved at distances larger than √ 8t = cl.
Finally, the evaluation of the coupling can be restricted to sectors of definite topological charge Q, cf. Eq. (4.6), leading to: In the remaining of this section we will determine the expectation value of n + n considering two different scenarios, a dilute gas of ordinary instantons with S 0 = 8π 2 , and one made of fractional instantons of action S 0 = 8π 2 /N .

E.1 Ordinary instantons
The basic asumption in the dilute gas approximation is that the gas is made of objects with positive (+1) and negative (−1) topological charge that are otherwise equally probable; the probality to create either of them per unit volume is given by R. Under this hypothesis, it is easy to determine the partition function restricted to the sector of topological charge Q, it is given by: and where I Q (x) stands for the modified Bessel function of the first kind. It is easy to evaluate n + n Q in this approximation, it is given by: We can also compute the topological susceptibility in this approximation. The partition function in the θ vacuum is given by: Z (θ ) = C n,n 1 n!n! (RV ) n+n e i(n−n)θ = Ce 2RV cos θ , (E.9) and the topological susceptibility can be computed from: 10) In the dilute gas approximation R is determined from the determinant of the fluctuation operator around one instanton. The result at one-loop order in infinite volume reads where n c is the number of collective coordinates, equal to 4N for ordinary SU (N ) instantons. The formula up to two loops including the prefactor has been computed in Ref. [81]. The resulting expression is infrared divergent, although on a finite volume it would be cutoff by the size of the box. Note however, that for the determination of λ dig Q we only need the relation x = 2RV = χ V , where one can use the measured topological susceptibility to determine the input parameter x and evaluate the coupling using eq. (E.8).

E.2 Fractional instantons
In the case of fractional instantons, the only difference stems from the fact that instantons contribute with charge ±1/N to the partition function. Note that the twist used in our simulations is in the class of the so-call orthogonal twists, satisfying n μνñμν /4 = 0 (mod N ) and therefore the total topological charge is still quantized in integer units. This fact has to be taken into account when formulating the dilute gas partition function for fractional instantons and therefore: We can now use that: I α (x) = n 1 n! (n + α + 1) leading to the following expression for the N = 3 topological susceptibility: The relation between x = 2RV and χ V is more involved in this case than for ordinary instantons, but nevertheless it can be numerically inverted to obtain x in terms of χ V .
In the dilute gas approximation, see i.e. [72,73]: where ρ is the characteristic fractional instanton size. The peculiarity of fractional instantons vs ordinary ones is that the collective coordinates only include the instanton location; in particular the size is not a moduli parameter and on a finite volume it is fixed in terms of the volume of the box. Therefore, we expect ρ ∝l in the small volume regime and ρ ∝ 1/ QC D for large volumes. When ρ ∝l, and setting μ = 1/(cl), we can write: 4 (λ(μ)) −2 e −8π 2 /λ(μ) μ=(cl) −1 , (E. 21) and determine x, up to an overall normalization constant, from the measurement of the TGF coupling.