Large-order NSPT for lattice gauge theories with fermions: the plaquette in massless QCD

Numerical Stochastic Perturbation Theory (NSPT) allows for perturbative computations in quantum field theory. We present an implementation of NSPT that yields results for high orders in the perturbative expansion of lattice gauge theories coupled to fermions. The zero-momentum mode is removed by imposing twisted boundary conditions; in turn, twisted boundary conditions require us to introduce a smell degree of freedom in order to include fermions in the fundamental representation. As a first application, we compute the critical mass of two flavours of Wilson fermions up to order $O(\beta^{-7})$ in a ${\mathrm{SU}}(3)$ gauge theory. We also implement, for the first time, staggered fermions in NSPT. The residual chiral symmetry of staggered fermions protects the theory from an additive mass renormalisation. We compute the perturbative expansion of the plaquette with two flavours of massless staggered fermions up to order $O(\beta^{-35})$ in a ${\mathrm{SU}}(3)$ gauge theory, and investigate the renormalon behaviour of such series. We are able to subtract the power divergence in the Operator Product Expansion (OPE) for the plaquette and estimate the gluon condensate in massless QCD. Our results confirm that NSPT provides a viable way to probe systematically the asymptotic behaviour of perturbative series in QCD and, eventually, gauge theories with fermions in higher representations.


Introduction
The success of perturbation theory in High Energy Physics (HEP) can hardly be denied. In particular, in asymptotically free theories, field correlators at short distances are reliably approximated by perturbative expansions in the running coupling at a large momentum scale. At the same time, even in these (lucky) cases, it is mandatory to have some control on nonperturbative effects, i.e. contributions that scale like powers of the QCD scale QCD . We will often refer to these as power corrections. A tool to take the latter into account was suggested back in the late seventies. This goes under the name of QCD sum rules, or Shifman-Vainshtein-Zakharov (SVZ) sum rules [1,2]. One of the authors defined the method as "an expansion of the correlation functions in the vacuum con-densates" [3]. These condensates are the vacuum expectation value of the operators that emerge in the Operator Product Expansion (OPE) for the relevant correlation function. In the OPE formalism the condensates are fundamental quantities, which are in principle supposed to parametrise power corrections in a universal way. By determining the value of a condensate in one context, one gains insight into different physical processes; this has in turn motivated several approaches to the determination of condensates. Having said all this, the sad news is that not all the condensates have actually the same status. In particular not all the condensates can be defined in a neat way, which ultimately means disentangled from perturbation theory. While this is the case for the chiral condensate, the same cannot be said for the gluon condensate, which is the one we will be concerned with in this work.
Based on a separation of scales, the OPE makes pretty clear what can/must be computed in perturbation theory, i.e. the Wilson coefficients. Still, this does not automatically imply that perturbative and nonperturbative contributions are separated in a clear-cut way. The key issue is that perturbative expansions in HEP are expected to be asymptotic ones on very general grounds. In particular, the series in asymptotically free theories are plagued by ambiguities which are due to so-called infrared renormalons [4,5]. From a technical point of view, renormalons show up as singularities which are encountered if one tries to Borel resum the perturbative series. All in all, there is a power-like ambiguity in any procedure one can devise in order to sum the series, and this ambiguity unavoidably reshuffles perturbative and nonperturbative contributions in the structure of the OPE. Being the Wilson coefficients affected by ambiguities that are power corrections, the general strategy is to reabsorb the latter in the definition of the condensates. This amounts to a prescription to give a precise meaning both to the perturbative series and to the condensates that appear in the OPE.
The idea of determining the gluon condensate from nonperturbative (Monte Carlo) measurements in lattice gauge theories dates back to the eighties and early nineties [6][7][8][9]. Based on symmetry grounds and dimensional counting, the two leading contributions in the OPE for the basic plaquette are given by the identity operator and the gluon condensate. Both operators appear multiplied by Wilson coefficients that can be computed in perturbation theory, and in particular the coefficient that multiplies the identity operator is simply the perturbative expansion of the plaquette. Other operators that appear in the OPE are of higher dimension, and their contributions are therefore suppressed by powers of a QCD . Subtracting from a nonperturbative (Monte Carlo) measurement of the plaquette the sum of the perturbative series, and repeating the procedure at different values of the coupling, the signature of asymptotic scaling, i.e. the signature of a quantity of (mass) dimension four, should become visible. With renormalons attracting more and more attention, it eventually became clear that such a procedure must be deeply affected by the ambiguities we discussed above, suggesting that a precise definition of the resummed perturbative expansion is necessary.
In the meantime Numerical Stochastic Perturbation Theory (NSPT) [10] was developed as a new tool for computing high orders in lattice perturbation theory. NSPT paved the way to the evaluation of many more terms in the perturbative expansion of the plaquette, and in turn made it at least conceivable that the behaviour of the series could be understood at the level of pinning down the correct order of magnitude of the ambiguity involved. Results of early investigations [11] were interesting: for the first time, it was clear that very high order contributions can be computed in perturbative series for lattice gauge theories. Unfortunately the pioneering NSPT studies of that time were far away from computing the series up to the orders at which the renormalon growth actually shows up in its full glory. With limited computing power available, a way out was sought in the form of a change of scheme (i.e. a scheme in which the renormalon behaviour is best recognised, possibly at lower orders than in the lattice scheme). Still, the numerical results were in the end puzzling as for consequences, since trying to sum the series from the information available even suggested the idea that an unexpected contribution from a dimension-2 operator was present [12]. Other attempts were made [13], but it eventually took roughly twenty years before the renormalon behaviour was actually captured [14][15][16][17], needless to say, via NSPT. 1 In SU(3) Yang-Mills theory the IR renormalon was indeed directly inspected, and the finite-size effects that are unavoidable on finite lattices assessed. The bottom line is that the victory is twofold. On one side, the renormalon growth is indeed proved to be present as conjectured (ironically, in a scheme -the lattice -which one would have regarded as the very worst to perform the computations). Given this, one has a prescription to sum the series and perform the subtraction (if sufficiently high orders are available, one can look for the inversion point in the series, where contributions start to grow and a minimum indetermination in summing the series can be attained).
The present work is a first attempt at performing the determination of the gluon condensate from the plaquette in full QCD, i.e. with fermionic contributions taken into account. The main focus here is in developing the NSPT technology, and present a first set of results, which allow a definition of the gluon condensate. In particular for the first exploration, we use existing Monte Carlo simulations for the plaquette in full QCD, as detailed below. Having ascertained that the pro-cedure is viable, a precise determination of the condensate in full QCD will require a dedicated Monte Carlo simulation, with a careful choice of the fermionic action. On top of being interesting per se, the methodology presented here opens the way to other applications, in which different colour groups and different matter contents can be investigated. The final goal would be to inspect whether in a theory that has an IR fixed point, the renormalon growth is tamed, as one would expect in theories where the condensates vanish. We defer these questions to future investigations, hoping to gain extra insight into the task of identifying the boundaries of the conformal window.
The paper is organised as follows. In Sect. 2 we review briefly how NSPT can be applied to lattice gauge theories. In Sect. 3 twisted boundary conditions for fermions in the fundamental representation are introduced. In Sect. 4 we discuss how to take into account fermions with smell in NSPT. We present our results for the expansion of the critical mass of Wilson fermions in Sect. 5, and for the expansion of the plaquette with staggered fermions in Sect. 6. In Sect. 7 we investigate the asymptotic behaviour of the expansion of the plaquette and extract the gluon condensate in massless QCD. In Sect. 8 we draw our conclusions and present some possible future steps.

Lattice gauge theories in NSPT
Let us here summarise the main steps in defining NSPT for lattice gauge theories. Rather than trying to give a comprehensive review of the method, we aim here to introduce a consistent notation that will allow us to discuss the new developments in the rest of the paper. For a more detailed discussion of the NSPT formulation, the interested reader can consult e.g. Ref. [18], whose notation we shall try to follow consistently. 2 In particular, we assume to work with a hypercubic lattice with volume L 4 = a 4 N 4 and assume the lattice spacing a to be 1, unless where stated otherwise. We use x, y, z for position indices, μ, ν, ρ = 1, . . . , 4 for Lorentz indices and α, β, γ = 1, . . . , 4 for Dirac indices.
The original formulation of NSPT is based on the Stochastic Quantization formulation of lattice field theories, in the case at hand lattice gauge theories. For the purposes of this study, we focus on gauge theories that are defined by the Euclidean Wilson action for the gauge group SU(N c ): where U is the product of the link variables, denoted U μ (x), around the 1 × 1 plaquette , and the sum extends over all the plaquettes in the lattice. Introducing a stochastic time t, a field U μ (x; t) can be defined that satisfies the Langevin equation As detailed in Appendix A, we have denoted by ∇ xμ the left derivative in the group; η is a stochastic variable defined in the algebra of the group, where T a are the generators of the group, and η a μ (x; t) are Gaussian variables such that The key point of Stochastic Quantization is that the larget distribution of observables built from the solution of the Langevin equation above corresponds to the distribution that defines the path integral of the quantum theory [19,20]: In order to develop NSPT, the dynamical variables U μ (x; t) can be expanded in powers of the coupling constant g, which is given in the lattice formulation by β −1/2 : Solving the Langevin equation, Eq. (2), order by order in β −1/2 yields a system of coupled equations for the perturbative components of the link variables U (k) μ (x; t). Expanding the solution of Langevin equation in powers of the coupling is a standard approach to proving the equivalence of stochastic and canonical quantisation, i.e. Eq. (5) [21], and was the starting point for stochastic perturbation theory: with this respect NSPT is just the numerical implementation of the latter on a computer. The idea of studying the convergence properties of a stochastic process order by order after an expansion in the coupling is actually quite general. In this spirit different NSPT schemes can be set up, also based on stochastic differential equations different from Langevin [22,23].
Euler integrator Discretising the stochastic time in steps of size allows a numerical integration of the Langevin equation, where the force driving the evolution is (8) and the operator g projects on the algebra (see Appendix A). Note that Eq. (8) does not lend itself to a perturbative solution in powers of β −1/2 , since there is a mismatch between the deterministic drift term, which starts at order β 1/2 , and the stochastic noise, which is of order β 0 . This is easily resolved by rescaling the integration step by a factor of β, so that both contributions start at order β −1/2 . Denoting the new time step τ = β, the force term becomes Expanding F in powers of β −1/2 , leads to a system of coupled equations for the evolution of the coefficients of the perturbative expansion of U . Omitting Lorentz and position indices, we get where η only contributes to the F (1) term.

Stochastic gauge fixing
The zero modes of the gauge action do not generate a deterministic drift term in the Langevin equation, and therefore their evolution in stochastic time is entirely driven by the stochastic noise, which gives rise to diverging fluctuations. This phenomenon is well known since the early days of NSPT, see e.g. Ref. [24], and is cured by the so-called stochastic gauge fixing procedure [25] applied to the theory formulated on the lattice. The procedure implemented in this work alternates an integration step as described above with a gauge transformation: where the field w(x) is defined in the algebra of the group, α is a free parameter, which we choose equal to 0.1 and ∇ * μ is the backward derivative in direction μ. Note that there is nothing compelling in the choice for w(x). In this work we make the same choice as in Ref. [24], which is slightly different from the one adopted in Ref. [18]: the corresponding gauge transformation does not lead, if iterated, to the Landau gauge. In NSPT the gauge transformation is expanded in powers of the coupling, and the transformation in Eq. (12) is implemented order by order in perturbation theory. The combined step for the integrator adopted in this work can be summarised as where all the terms are expanded in powers of β −1/2 , and the perturbative components are updated.

Runge-Kutta integrator
Higher order integrators, in particular Runge-Kutta schemes, have been used for the lattice version of the Langevin equation since the early days [20]. A new, very effective second-order integration scheme for NSPT in lattice gauge theories has been introduced in Ref. [15]. While we have tested Runge-Kutta schemes ourselves for pure gauge NSPT simulations, in this work we adhere to the simpler Euler scheme: when making use of the (standard) stochastic evaluation of the fermionic equations of motion (see Sect. 4), Runge-Kutta schemes are actually more demanding (extra terms are needed [26,27]).

Twisted boundary conditions and smell
When a theory is defined in finite volume, the fields can be required to satisfy any boundary conditions that are compatible with the symmetries of the action. We adopt twisted boundary conditions (TBC) [28] in order to remove the zeromode of the gauge field, and have an unambiguous perturbative expansion, which is not plagued by toron vacua [29]. The gauge fields undergo a constant gauge transformation when translated by a multiple of the lattice size; therefore twisted boundary conditions in directionν are where μ ∈ SU(N c ) are a set of constant matrices satisfying Fermions in the adjoint representation can be introduced in a straightforward manner; the boundary conditions with the fermionic field in the matrix representation read The inclusion of fermions in the fundamental representation is not straightforward; indeed, the gauge transformation for the fermions when translated by a multiple of the lattice size reads leading to an ambiguous definition of ψ(x + Lμ + Lν). An idea to overcome this problem, proposed in Ref. [30] and implemented e.g. in Ref. [31], is to introduce a new quantum number so that fermions exist in different copies, or smells, which transform into each other according to the antifundamental representation of SU(N c ). The theory has a new global symmetry, but physical observables are singlets under the smell group. Thus, configurations related by a smell transformations are equivalent, and in finite volume we are free to substitute Eq. (19) with where ν ∈ SU(N c ). It is useful to think of the fermion field as a matrix in colour-smell space. If the transformation matrices in smell space satisfy the same relations as in Eq. (17) (in particular we choose them to be equal to the s), then twisted boundary conditions are well-defined.
It is worth pointing out that, through a change of variable in the path integral [32,33], twisted boundary conditions could be equivalently implemented by multiplying particular sets of plaquettes in the action by suitable elements of Z N c and considering the fields to be periodic. This change of variable works only in the pure gauge or fermions in the adjoint representation cases. Thus, the explicit transformation of Eq. (20) is required when fermions in the fundamental representation with smell are considered.

Fermions in NSPT
is the action of a single fermion, then dynamical fermions in NSPT can be included thanks to a new term in the drift, as shown in Refs. [20,34]: the determinant arising from N f degenerate fermions can be rewritten as and can be taken into account by adding −N f Tr ln M to the gauge action. From the Lie derivative of the additional term and recalling that a rescaled time step τ = /β is used in the Euler update, we obtain the new contribution to be added to the pure gauge drift. It is important to note that the coefficient of i T a is purely real because the Wilson oper-ator is γ 5 -Hermitian and the staggered operator is antihermitian: this is consistent with the drift being an element of the algebra. The trace can be evaluated stochastically: Eq. (22) is replaced by thanks to the introduction of a new complex Gaussian noise ξ satisfying 3 ξ * (y) βir ξ(z) γ js = δ yz δ βγ δ i j δ rs .
The real part must be enforced, otherwise the dynamics would lead the links out of the group since the drift would be guaranteed to be in the algebra only on average. In NSPT, the Dirac operator inherits a formal perturbative expansion from the links, M = ∞ n=0 β −n M (n) , so the inverse ψ = M −1 ξ can be computed efficiently from the knowledge of the inverse free operator via the recursive formula The inverse of the free operator is conveniently applied in Fourier space. If fermions have smell, then the rescaling N f → N f /N c is required in order to have N f flavours in the infinite-volume limit. In other words, this is the same as considering the N c th root of the determinant of the fermion operator. In principle such rooted determinant could come from a nonlocal action, because twisted boundary conditions break the invariance under smell transformations. Nevertheless, this rooting procedure is sound since we know in advance that in the infinite-volume limit all the dependence on boundary conditions will be lost and the determinant will factorise as the fermion determinant of a single smell times the identity in smell space. It is also possible to show with arguments similar to those presented in Ref. [35] that, if the theory without smell is renormalisable, this operation leads to a perturbatively renormalisable theory as well. Below we describe in detail Wilson and staggered fermions in the fundamental representation, so we explicitly rescale N f → N f /N c . It is also important to remember that the fermion field, seen as a matrix in colour-smell space, is not required to be traceless, thus its Fourier zero-mode does not vanish: we require antiperiodic boundary conditions in time direction not to hit the pole of the free propagator in the massless case. We avoid twisted boundary conditions in time direction because in the massless case it might happen for the free fermion propagator to develop a pole at some particular momenta.

Wilson fermions
The Wilson Dirac operator and its Lie derivative are M yβir,zγ js = (m + 4)δ rs δ yz δ βγ δ i j where the non-diagonal term has been expressed through We must give a perturbative structure to the mass m = ∞ n=0 β −n m (n) to account for an additive mass renormalisation, see Sect. 5. The stochastic evaluation of the trace leads to where ϕ (μ) = D(μ)ψ,φ (μ) = γ 5 D(μ)γ 5 ξ and the fermion fields have been represented as matrices in colour-smell space. After taking the real part, the fermion drift can be finally written as In Appendix B the actual implementation of the fermion drift is described (only one of the two terms in Eq. (29) is actually needed). With the Fourier transform described in Appendix C, the inverse free Wilson operator with twisted boundary conditions is diagonal in momentum space and can be expressed as

Staggered fermions
We implemented for the first time staggered fermions in NSPT. The staggered field has no Dirac structure and describes four physical fermions in the continuum limit. Therefore, we rescale N f → N f /4 and the staggered operator is understood to be rooted when the number of flavour is not a multiple of four. The staggered Dirac operator and its Lie derivative are where the non-diagonal term has been expressed through and x ν is the staggered phase. The stochastic evaluation of the trace is analogous to the Wilson fermion case and Eq. (28) becomes with ϕ (μ) = D(μ)ψ andφ (μ) = −D(μ)ξ , leading to the final expression Again, the actual implementation of the staggered drift is shown in Appendix B. With the Fourier transform described in Appendix C, the inverse free staggered operator with twisted boundary conditions is found to be where1 = 0, μ + 1 =μ+μ andδ is the periodic Kronecker delta, with support in 0 mod 2π . The propagator is not diagonal in momentum space because the action depends explicitly on the position through α μ (x), but it is simple enough to avoid a complete matrix multiplication over all the degrees of freedom. If we aim to compute M (0) −1 v for some field v in momentum space, it is useful to represent v( p ) p ⊥ as matrices N c × N c with indicesñ 1 ,ñ 2 defined at each p site (n 1 , n 2 , n 3 , n 4 ) (see again Appendix C). Then the nondiagonal terms become diagonal when shifting iteratively v by L/2 in the p space. Incidentally, we must consider L to be even so that at the same time L/2 is well defined and (in the massless case) no spurious pole is hit when Eq. (35) is evaluated in finite volume: this stems from the fact that the staggered action is only invariant under translation of two lattice spacings, therefore twisted boundary conditions would be inconsistent for L odd.

The critical mass of Wilson fermions
The inverse of the Wilson fermion propagator in momentum space can be expressed as is the self energy. In this section the lattice spacing a is written explicitly. Wilson fermions are not equipped with chiral symmetry when the bare mass m vanishes: the self energy at zero momentum is affected by a power divergence a −1 , which has to be cured by an additive renormalisation. In an on-shell renormalisation scheme, the critical value of the bare mass, m c , for which the lattice theory describes massless fermions, is given by the solution of As observed in Ref. [36], this prescription matches the one obtained by requiring the chiral Ward identity to hold in the continuum limit. Expanding Eq. (37) defines the critical mass order by order in perturbation theory. The perturbative expansion of the inverse propagator is where we have indicated explicitly the dependence of the coefficients on the bare mass am. The functions (n) (ap, am) are matrices in Dirac space; since we are interested in the small momentum region and (n) (0, am) is proportional to the identity, we consider (n) (ap, am) as scalar functions: when ap = 0 a projection onto the identity is understood.
Plugging the perturbative expansion of the critical mass where the dependence of γ (n) on m (n) c has been made explicit and c . Therefore, the renormalisation condition in Eq. (37) becomes order by order For illustration, we can compute the recursive solution of Eq. (37) at first-and second-order in the expansion in powers of β −1 , which yields Both results are familiar from analytical calculations of the critical mass. The first equation encodes the fact that the mass counterterm at first order in perturbation theory is given by the one-loop diagrams computed at zero bare mass. The second equation states that the second-order correction is given by summing two-loop diagrams evaluated at vanishing bare mass, and one-loop diagrams with the insertion of the O β −1 counterterm, see e.g. Ref. [37].
It should also be noted that, when working in finite volume, momenta are quantised. Unless periodic boundary conditions are used, p = 0 is not an allowed value for the momentum of the states in a box. Therefore, condition (37) can only be imposed after extrapolating the value of to vanishing momentum. The detailed implementation is discussed below in Sect. 5.1.
Critical masses have been computed analytically up to two loops [37,38], and in NSPT at three and four loops [39,40]. High-order perturbation theory with massless Wilson fermions requires the tuning of the critical mass at the same order in β −1 , and it is possible to determine this renormalisation using NSPT. Let us illustrate the strategy in detail. We begin by collecting configurations for different time steps τ of the stochastic process; for each configuration the gauge is fixed to the Landau gauge [41,42]. The propagator at momentum p is computed by applying the inverse Dirac operator to a point source in momentum space, For each simulation at a given value of τ , the error bars are computed as detailed in Appendix D. The propagator with periodic boundary conditions is a (diagonal) matrix in colour and momentum space and has a Dirac structure; it is important to stress again that with TBC there is not a colour structure any more and the momentum has a finer quantisation. The average over all the configurations gives the Monte Carlo estimate of S( p). We can now extrapolate the stochastic time step to zero and invert the propagator to obtain S( p) −1 . Finally, the inverse propagator is projected onto the identity in Dirac space. All these operations are performed order by order in perturbation theory keeping in mind that, after the measure of the propagator, all perturbative orders β −k/2 with an odd k are discarded, since the expansion in powers of β −1/2 is an artefact of NSPT. The errors can be estimated by bootstrapping the whole procedure. The legacy of this process is the measure of the functions γ (n) (ap), as it is clear from Eq. (40). The renormalisation condition in Eq. (41) must then be imposed: this can be done iteratively one order after the other. When all the coefficients up to some m (n) c are included in the simulation, all the γ functions up to γ (n) (ap) extrapolate to zero; on the other hand, from . In order to move on and compute the following coefficient of the critical mass, a new set of configurations where m (n+1) c is taken into account must be generated.
The procedure we described is well defined and even theoretically clean, since it enlightens the status of our m c as a perturbative additive renormalisation: once it is plugged in at a given order, the renormalised mass turns out to be zero at the prescribed order. On the other side, it is not at all the only possible procedure. The prescription of the authors of Ref. [23] is to expand the solution of the stochastic process both in the coupling and in the mass counterterm. This is in the same spirit of Ref. [43]: the solution of the stochastic process can be expanded in more than one parameter and once a precise power counting is in place, the resulting hierarchy of equations can be exactly truncated at any given order. There are pros and contras for both approaches, i.e. the one we followed and the double expansion. The latter can provide a better handle on estimating errors due to the critical mass value; on the other side, it is expected to be numerical more demanding. All in all, we did not push Wilson fermions to very high orders: moving to the staggered formulation was by far the most natural option for the purpose of this work.

Zero-momentum extrapolation and valence twist
Since in finite volume it is possible to measure (ap) only for discretised non-zero momenta, the data need to be extrapo-lated to zero momentum using a suitable functional form. The strategy adopted in the literature -see for example Eqs. (13) and (14) in Ref. [40] -is based on expanding the quantities of interest in powers of ap. In the infinite-volume limit, such an expansion leads to a hypercubic symmetric Taylor expansion composed of invariants in ap, logarithms of ap and ratios of invariants; an explicit one-loop computation to order a 2 is shown e.g. in Eq. (24) of Ref. [44]. The ratios and the logarithms arise because we are expanding a nonanalytic function of the lattice spacing: infrared divergences appear when expanding the integrands in ap. On the other hand, working consistently in finite volume does not cause any infrared divergence: expressions for γ (n) (ap) will be just sums of ratios of trigonometric functions, which we can expand in ap obtaining simply a combination of polynomial lattice invariants. 4 Still, this is not enough for a reliable extrapolation to vanishing momenta. In order to understand better the range of momenta that allow a reliable extrapolation, we computed γ (1) (ap) in twisted lattice perturbation theory (see Appendix E). As a cross-check of our calculation we verified that γ (1) (0) is gauge-invariant (this result must be true at all orders because of the gauge-invariance of the pole mass [45]). It can be seen from the analytic expansion of γ (1) (ap) that even the lowest momentum allowed on our finite-size lattices, ap 1,2,3 = 0, ap 4 = π/L, is far from the convergence region of this series. This happens even for reasonably big lattices, L 32. In order to increase the range of available momenta, we use θ -boundary conditions [46] for the valence fermions, thereby reaching momenta p 4 = θ/L which are within the convergence radius of the ap-expansion. The hypercubic series becomes just a polynomial in (ap 4 ) 2 by setting all the other components to zero. The agreement between data and the analytic finitevolume calculations can be seen in Fig. 1. It is worthwhile to emphasise that measuring such low momenta requires a careful analysis of the thermalisation. At the lowest order we can check directly when the measures agree with the theoretical predictions. At higher orders, it is necessary to wait until the statistical average has clearly stabilised, as shown in Fig. 2. This kind of analysis is computationally intensive: in the case at hand, we performed up to 5 · 10 6 lattice sweeps, saving one propagator every 10 3 sweeps. The first 2 · 10 3 configurations have been discarded in the analysis. is used for fitting. Most analytic finite-volume predictions have been drawn as lines to help the eye in the comparison. The difference with the prediction in the right panel is to be ascribed to the fact that we are able to resolve finite volume effects We determined the first 7 coefficients of the critical mass for N c = 3 and N f = 2 on a 16 4 lattice with twisted boundary conditions on a plane. The twist matrices are corresponding to z 12 = exp i 2π 3 . Configurations are collected at three different time steps, τ = 0.005, 0.008, 0.01. Because the volume and the number of colours are large compared to the former test in Fig. 1, it is computationally too expensive to replicate the same statistics at all orders: we settled for 5 · 10 5 sweeps at the smallest τ , measuring the propagator every r = 10 3 sweeps. At larger time steps, we rescale these numbers to keep the product r · τ constant. The propagator is measured at the smallest available momentum, which has θ/L in the time component and vanishes elsewhere; we choose three different values for the phase of the valence twist, θ = π/2, 2π/3, 4π/5. Extrapolations to zero momentum are performed using a linear fit in (ap) 2 . The analysis is performed on different subsets of the data 5 to estimate systematic errors. The total error is the sum in quadrature of half the spread around the central value among the different fits and the largest error from the fits.
The procedure described in Sect. 5.1, even though welldefined, is found to be numerically unstable at high orders. The number of propagators required to reach a clear plateau, like the ones shown in Fig. 2, is beyond what it can be reasonably collected with the current NSPT implementations. Therefore, we decided to proceed with a smaller statistics and to add a new systematic uncertainty for the extrapolated coefficients, as explained below. It has to be emphasised that once a coefficient of the critical mass is determined, only the central value is used as input for the following runs: even if we could collect enough statistics and manage to reduce the error, that is not included in the simulations. This makes c . Although γ (1) (ap) does not extrapolate to zero, the extrapolation of γ (4) (ap) is compatible with the value known from Ref. [40]. Notation as in Fig. 1   Fig. 4 Determination of the coefficient m (8) c . The errors overshadow the value of the critical mass, which is compatible with zero. Notation as in Fig. 1 the impact of the uncertainty of m (n) c on m (n+1) c and higher hard to assess; also, performing simulations for several values of each coefficient is not feasible. To be conservative, we adopted the following strategy. Once a critical mass m (n) c is determined and put in the next-order simulation, the corresponding γ (n) (ap) should extrapolate to zero. If it extrapolates to n , we take | n /m (n) c | as an estimate of the relative systematic error to be added in quadrature to the determination of all the higher-order critical masses.
Despite these instabilities, the lower-order results are close to the known coefficients (keeping in mind that we might resolve finite-volume effects), as it can be seen for example in Fig. 3. We stopped the procedure at m (8) c , when the errors started dominating over the central value of the coefficient, see Fig. 4. Our results are summarised in Table 1.

Perturbative expansion of the plaquette
Following Ref. [16], we define the average plaquette so that the value of P ranges between 0, when all link variables are equal to the identity, and 1. The plaquette expectation value has the perturbative expansion the coefficients p n are obtained from the Langevin process.

Numerical instabilities
The study of the NSPT hierarchy of stochastic processes is not trivial. While there are general results for the convergence of the generic correlation function of a finite number of perturbative components of the fields [18,47], the study of variances is more involved, and many results can only come from direct inspection of the outcome of numerical simulations. In particular, one should keep in mind that in the context of (any formulation of) NSPT, variances are not an intrinsic property of the theory under study; in other words, they are not obtained as field correlators of the underlying theory. Big fluctuations and correspondingly huge variances were observed at (terrifically) high orders in toy models [47]: signals are plagued by several spikes and it is found by inspection that a fluctuation at a given order is reflected and amplified at higher orders. All in all, variances increase with the perturbative order (not surprisingly, given the recursive nature of the equations of motion). Moving to more realistic theories, a robust rule of thumb is that, as expected on general grounds, the larger the number of degrees of freedom, the less severe the problems with fluctuations are. In particular, we have not yet found (nor has anyone else reported) big problems with fluctuations in the computation of high orders in pure Yang-Mills theory. We now found that the introduction of fermions indeed causes instabilities at orders as high as the ones we are considering in this work. Once again, this effect can be tamed by working on increasingly large volumes. Once a fluctuation takes place, the restoring force would eventually take the signal back around its average value but in practice this mechanism is not always effective. At high orders the instabilities can be so frequent and large that the signal is actually lost, and the average value of the plaquette becomes negligible compared to its standard deviation, as it is illustrated in Fig. 5. The order at which the signal is lost is pushed to higher values by increasing the volume, but eventually uncontrolled fluctuations will dominate. Moreover, we find that spikes tend to happen more frequently at smaller τ . Roughly speaking, this does not come as a surprise, since at smaller time steps one has to live with a larger number of sweeps, thereby increasing the chances of generating large fluctuations when computing the force fields. In Table 2 the orders available at each volume and time step are shown in detail.

Determination of the p n
The lowest coefficients have already been computed analytically. In particular, in twisted lattice perturbation theory we have that is volume independent [48]. The infinite-volume value of p 1 can be obtained adding to the pure gauge contribution [49], the contribution due to staggered fermions [50], For the specific case N c = 3, N f = 2, we find p 1 = 1.10312 (7). We also computed the fermion contribution to   (6) for N f = 2. Trying to extract p 0 and p 1 from our data at L = 48, we realise that even τ 2 effects in the extrapolation must be considered because of the very high precision of the measurements. For these two coefficients, a dedicated study at has been performed, which required new simulations at time steps τ = 0.004 and τ = 0.0065; the agreement with the analytic calculations is found to be excellent, see Fig. 6. Therefore, p 0 and p 1 are set to their infinite-volume values and excluded from the analysis of the numerical simulations. The remaining orders are obtained from NSPT. The value p n,τ for the plaquette at order n and time step τ is computed from the average of the fields generated by the stochastic process, after discarding a number of thermalisation steps. The moving averages result to be stable, as can be seen in the two examples of Fig. 7. In order to exploit all the available data, the thermalisation is set differently at different orders. The covariance Cov(n, m) τ between p n,τ and p m,τ is computed taking into account autocorrelations and cross-correlations, as explained in detail in Appendix D. Clearly there is no correlation between different τ . In order to estimate the covariance when two orders have different thermalisations, we take into account only the largest set of common values where both are thermalised. This pairwise estimation of the covariance matrix does not guarantee positive definiteness, therefore we rely on Higham's algorithm, where the coefficients a n are the slopes of the combined linear fits. The interesting fit results are the values of the extrapolated plaquettes p n and their covariance matrix Cov(n, m).
In general, because of the available statistics and the intrinsic fluctuations of the observable, the lower-order values are measured more accurately compared to the higher-order ones; the same holds for the estimate of the entries the covariance matrix. Since, in principle, the plaquette at a certain order could be extracted without any knowledge about its higher-order values, we can get the best estimate for a p n by implementing the fit iteratively, increasing n max from 0 to the maximum available order. At each iteration, we determine the order with the minimum number of measures N min and rescale the entries of the covariance matrix so that there is a common normalisation (N = N min in Eq. (83)) for all the matrix elements. In this way, all the data are exploited for the determination of the covariance of the process, and the non-positive definiteness of the covariance of the averages arises only from the presence of autocorrelations and cross-correlations. Higham's algorithm is then applied to Cov(n, m) τ restricted to n max orders. At this stage, minimising the χ 2 allows us to extract p n max with Cov(n max , m) for m ≤ n max . The tolerance of Higham's algorithm is tuned so that the covariance matrix is able to represent our data, i.e. so that the reduced chi-squared is close to 1. The combined fit determines also the plaquettes at orders lower than n max , which are always checked and found to be in agreement, within errors, with their previous determination at smaller n max . An example of a correlation matrix extracted with this procedure is in Fig. 8, where clear structures of correlated and anticorrelated coefficients are visible. The results of the combined extrapolations are summarised in Table 3.

Gluon condensate
In this section we restore the lattice spacing a and follow the notation of Refs. [16,17]: the gluon condensate is defined as the vacuum expectation value of the operator where the coupling α is related to the Wilson action coupling by α = N c 2πβ and the beta function is with the scheme-independent coefficients It is useful to remember that, in the massless limit, O G is renormalisation group invariant and depends on the scheme only through the renormalisation condition used to define the composite operator. It is easy to relate the gluon condensate and the plaquette in the naive continuum limit:   In the interacting theory mixing with operators of lower or equal dimension occurs. For the case of the plaquette, the mixing with the identity needs to be considered, yielding (56) which shows explicitly the subtraction of the quartic power divergence. 7 As a consequence where P MC is the plaquette expectation value obtained from a nonperturbative Monte Carlo simulation. As such, P MC is expected to depend on the cut-off scale a, and QCD . In the limit a −1 QCD , Eq. (57) can be seen as an Operator Product Expansion (OPE) [1,2,53], which factorises the dependence on the small scale a. In this framework, 8 condensates like O G are process-independent parameters that encode the nonperturbative dynamics, while the Wilson coefficients are defined in perturbation theory, Note that both Z and C G depend only on the bare coupling β −1 , and do not depend on the renormalisation scale μ, as expected for both coefficients [55,56]. Nonperturbative contributions to Z , or C G , originating for example from instantons, would correspond to subleading terms in QCD . This procedure defines a renormalisation scheme to subtract power divergences: condensates are chosen to vanish in pertubation theory or, in other words, they are normal ordered in the perturbative vacuum. This definition matches the one that is natural in dimensional regularisation, where power divergences do not arise. Nevertheless, it is well known that such a definition of the condensates might lead to ambiguities, since the separation of scales in the OPE does not necessarily correspond to a separation between perturbative and nonperturbative physics (see the interesting discussions in Refs. [3,57]). For example, the fermion condensate in a massless theory is well-defined since, being the order parameter of 7 We mention that, in a theory with fermions, the operator O G must be combined with mψψ to give a renormalisation group invariant quantity; moreover mixing with the operators mψψ andψ(i / D − m)ψ should also be considered [51,52]. Clearly such complications are not present in the massless case and the operator iψ / Dψ can be neglected in the following discussions since it vanishes when the equation of motion are used. 8 It is useful to keep in mind that other definitions of the gluon condensate are possible, see e.g. Ref. [54]. chiral symmetry breaking, it must vanish in perturbation theory. The same cannot be said for the gluon condensate [58], and indeed the ambiguity in its definition is reflected in the divergence of the perturbative expansion of the plaquette. For this picture to be consistent, it must be possible to absorb in the definition of the condensate the ambiguity in resumming the perturbative series.
In the following, we are going to study the asymptotic behaviour of the coefficients p n determined in the previous section and discuss the implications for the definition of the gluon condensate in massless QCD.

Growth of the coefficients
From the analysis in Refs. [11,16], it is possible to predict the asymptotic behaviour of the ratio where the use of the Wilson action with N c = 3 is assumed. This relation can be derived under the hypothesis that the plaquette series has a fixed-sign factorial divergence and the corresponding singularity in the Borel plane is the source of an ambiguity that can be absorbed by redefining the condensate. It is not possible to go further in the 1/n expansion since the β 2 coefficient is scheme-dependent and it is not known for staggered fermions. In Figs. 9 and 10, the comparison between Eq. (59) and our data at different volumes is shown. How finite-volume effects influence the values of the coefficients p n has already been studied in the literature [16,59]. From a standard renormalon-based analysis, the value of the loop momenta that contribute the most to p n decreases exponentially with n. Since the finite size of the lattice provides a natural infrared cutoff, we expect finite-volume effects to be larger at larger perturbative orders. The dependence of p n on the lattice size N can be modelled with a finite-volume OPE, exploiting the separation of scales a −1 (N a) −1 : the leading correction is [16] where α((N a) −1 ) must be expressed in terms of the coupling β at the scale a −1 using Eq. (53). We do not attempt to take into account 1/N 4 effects, as our data do not allow to perform a reliable combined fit. Apparently no significant finite-volume effects are visible where they would be expected the most, i.e. at larger n. This is shown in the two examples of Fig. 11. A similar behaviour has been observed Fig. 9 Ratio p n /(np n−1 ) extracted from our data at L = 24, 28, 32, 48. In order to be visible, points referring to different volumes are placed side by side. The leading order (LO) prediction refers to the n → ∞ limit, while the next-to-leading order (NLO) one includes the first 1/n correction Fig. 10 Same as Fig. 9, but the region at large n is enlarged in Ref. [16], where the data points computed on comparable volumes show little dependence on the lattice size. In that study, a detailed analysis with a large number of volumes was needed in order to be able to fit the finite-volume corrections. The overall effect is found to be an increase of the ratio p n /(np n−1 ), see e.g. Fig. 6 in Ref. [16]. In our case, data in finite volume do cross the theoretical expectation; still, considering the spread between points at different volumes in Fig. 10 as a source of systematic error, we could consider our measurements to be compatible with the asymptotic behaviour of Eq. (59). We also ascertain the existence of an inversion point when resumming the perturbative series, as explained in Sect. 7.3. Despite this encouraging behaviour, any definite conclusion about the existence of the expected renormalon can only be drawn after performing an appropriate infinite-volume study. We emphasise that in this work the discrepancies in the determination of the p n from different volumes must be interpreted as part of our systematic uncertainty, being this an exploratory study. A precise assessment of the finite-volume effects will be sought for a precise determination of the gluon condensate; we are currently planning a set of dedicated simulations in the near future to settle this issue.

Monte Carlo plaquette
Nonperturbative values for the SU(3) plaquette with N f = 2 (rooted) staggered fermions can be found in Ref. [50], where data are collected from Refs. [60,61]. For each value of the  The grey points are available from Ref. [50] but are excluded because of our fit procedure. In most cases the error bar is smaller than the symbol. The orders of the polynomials used in the fits are in Table 4 bare coupling, the physical scale is provided via the Sommer parameter r 0 [62]. The data are given for several values of the fermion bare mass, and need to be extrapolated to the chiral limit for our purposes. A reasonable assumption (for example adopted and verified also in Ref. [63] for the ratio r 0 /a) is that the plaquette and the ratio r 0 /a have a polynomial behaviour at small masses. We performed fits with linear to cubic polynomials and varied the fit ranges to exclude points at larger values of the masses, but in many cases the fits did not return a satisfactory description of the data with sensible values of χ 2 /dof. Because we are using results from past simulations, it is difficult to track accurately the systematic errors in the data. For this reason, we decided to choose the fit with smaller χ 2 /dof among those we tried and if χ 2 /dof > 1 the errors in the data were rescaled by a common factor in order to have a reduced chi-squared equal to 1. The fits resulting from this approach are shown in Fig. 12; the extrapolated values for plaquettes and r 0 /a are in Table 4. Another approach consists in considering the average between the largest and smallest extrapolated values among all the different fits we tried (without rescaled errors and with reduced chi-squared smaller than some reasonable threshold) and assigning an error equal to the sum in quadrature between the largest error from the fits and half the difference between the largest and smallest extrapolated values. In this way we obtain results compatible (both for central values and errors) with those in Table 4, confirming that the chiral extrapolation is sound and the error bars conservative enough. Note that in this paper we do not aim at a precise determination of the condensate, and therefore we can be satisfied with an inflated error on the Monte Carlo data points.

Determination of the minimal term
The perturbative contribution to the plaquette can be defined by the sum of the series up to the minimal term. The determination of the minimal term, and the summation of the series   are performed separately for each volume. We choose the prescription adopted in Ref. [17], i.e. we define the minimal term to be the valuen that minimises the product p n β −(n+1) and resum the series, Our results for all combinations of L and β are summarised in Table 5. The ordern at which the series starts to diverge depends only on the central value of the coefficients p n and not on their errors: in order to check that the inversion point determined by our procedure is stable, we bootstrapped the procedure by generating an ensemble of sets of coefficients { p n }. For each set, the coefficients p n are drawn from a Gaussian probability, whose mean and covariance are taken from the fit procedure described in Sect. 6. We then determinen for each of these sets. The inversion point turns out to be stable, as shown in Fig. 13 for a the case L = 48, and β = 5.3. This particular case is shown for illustration purposes, and the same features are seen in all other combinations of L and β.
The gluon condensate is then determined from with The coefficient β 2 is not universal, and is actually unknown for the discretisation used in this work. Not knowing β 2 prevents us from going further in the expansion of C G ; since the correction due to the Wilson coefficient falls between 5% and 6% for the values of β considered, a 6% systematic uncertainty is added in quadrature after the subtraction. The result of the subtraction is shown in the left panel of Fig. 14, for the largest volume. Since only a few values of β is available, it is hard to assess unambiguously the presence of a plateau. We decided to discard from the analysis the two values of the coupling corresponding to the coarser lattices, and define our best estimate of the condensate as the weighted average of the values obtained at the remaining βs. Our final results are summarised in the first column of Table 6. Table 6 Determination of the gluon condensate at different volumes. The determination labelled with 1 is obtained from the weighted average of the values at the three largest values of β. The determinations labelled with 2 and 3 are obtained by studying the scaling of a 4 O G with a 4 , as in the right panel of Fig. 14; they correspond respectively to the fit without and with a 6 correction (see text for the details) In order to put the choice of fit range on more solid ground, we studied the scaling of a 4 O G as a function of a 4 , as shown in Fig. 14. The slope of a linear fit of the three finest lattice spacings should give a determination of the condensate compatible with the value extracted from the weighted average. The spread between these two determinations and among the different volumes gives an idea of the magnitude of the systematic uncertainties involved. We also tried to include in the analysis all the available values of β and add a a 6 correction, in the attempt to model the deviations at large values of the coupling; this procedure gives again consistent results (despite a larger χ 2 ).
Truncating the sum up to the minimal term is one of the possible prescriptions to define the sum of a divergent series. The intrinsic ambiguity associated to S P (β) can be defined as the imaginary part of the Borel integral, which at leading order in 1/n is √ πn/2 pn β −n−1 [5]. In Table 7, the ambiguity associated to the gluon condensate is summarised. 9 9 Our definition of the ambiguity differs from the one in Ref. [16] by a factor √ π/2.

Conclusions
We used NSPT to perform for the first time large-order computations in lattice gauge theories coupled to massless fermions. We adopted twisted boundary conditions for the gauge fields to remove the zero-momentum mode. Since our fermions are in the fundamental representation, we consistently provided them with a smell degree of freedom. Both Wilson and (for the first time in NSPT) staggered fermions have been implemented. While for the former we performed an exploratory study of the critical mass up to order O(β −7 ), the latter are ultimately the best choice to reach very high orders, due to their residual chiral symmetry that bypasses the need of an additive mass renormalisation. Numerical instabilities were noticed in the study of simple models in NSPT since the early days of the method, but gauge theories have always been reported to stay on a safe side in this respect, even at orders as high as the ones we investigated in this work. With fermions in place, we now found that numerical instabilities arise for lattice gauge theories at high orders. While we plan to investigate the causes and develop a solution to this, the problem did not prevent us to reach order O(β −35 ) in the expansion of the basic plaquette for N c = 3 and N f = 2.
The plaquette has been for a long time the stage for the determination of the gluon condensate, to which is connected in the continuum limit. The perturbative expansion of the plaquette, which corresponds to the power divergent contribution associated to the identity operator in the relevant OPE, must be subtracted from nonperturbative Monte Carlo lattice computations. This long-standing and tough problem was eventually solved a few years ago in pure gauge [16,17], thanks to NSPT. Equipped with our high-orders expansions, we tackled once again the problem in the lattice regularisation of full QCD. We computed the perturbative expansion of the plaquette, and subtracted it from Monte Carlo measurements. In this context, NSPT is crucial: it is actually the only tool enabling this procedure, which asks for having the asymptotic behaviour of such series under control. This happens since the perturbative expansion of the plaquette is expected to be plagued by renormalon ambiguities. Under the assumption of considering finite-volume effects as a source of systematic errors, the observed growth of the coefficients in the expansion could be compatible with the leading IR renormalon; nevertheless, the large uncertainties and the lack of a study of finite-volume effects prevent us from drawing any definite conclusion. The IR renormalon forces to absorb the ambiguities attached to the perturbative series into the definition of the condensate itself. All in all, this implies that we needed a prescription to perform the computation. The one we chose amounts to summing the perturbative series up to its minimal term (which means computing the series up to orders that only NSPT can aim at). We regard this project as a first exploratory study. We could confirm both that the IR renormalon can be directly inspected, and that the series can be computed up to orders where the inversion point beyond which the expansion starts to diverge (at values of the coupling which are the typical ones in lattice simulations) is clearly visible. We performed our simulations at different lattice extents, in order to have a first estimate of finite-size effects (again, in both the study of renormalon behaviour and in the truncation of the series). This is the point which has to be better investigated in a following study. At the moment, finite-size effects are still to be considered as a systematic source of errors in our procedure.
On top of the follow-ups we have already discussed, we plan to extend our study to different number of colours, number of flavours and fermionic representations. It would be of the utmost importance to assess the high-order behaviour of perturbative coefficients in gauge theories different from QCD, to probe regions in the space of theories in which a (quasi-)conformal window can be present. This could be a powerful, alternative method to look for candidate theories for physics beyond the Standard Model.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

A Group theory conventions
The conventions used for group theoretical manipulations are summarised here. We consider the gauge group SU(N c ).
The generators of the group are denoted by T a ; the indices a, b, c = 1, . . . , N 2 c − 1 are assumed to be indices in the adjoint representation. The generators are defined to be Hermitian, and satisfy the commutation relations where f abc are the group structure constants. The normalisation of the generators is chosen to be such that The left derivative on the group is defined as ∇ xμ = a T a ∇ a U μ (x) , where the Lie derivative is given by We define an operator, g , that projects on the algebra g of the group: The indices i, j = 1, . . . , N c will be used as indices in the fundamental representation, r, s = 1, . . . , N c as indices in the antifundamental representation.

B Optimisation of the fermion drift
A useful optimisation consists in improving on Eqs. (29) and (33) so that it becomes numerically cheaper to evaluate the fermion drift. Considering for example Wilson fermions, we notice that it is possible to simplify the trace whereTr is tracing all indices but the position one, and we used the fact that the inverse Wilson operator is γ 5 -Hermitian. For staggered fermions the simplification is analogous because the inverse staggered operator is antihermitian. The step must be done before the stochastic evaluation of the trace: once the random sources are introduced, cyclic invariance gets broken and will be restored only on average. Using Eq. (69) as a starting point, we obtain a drift which is already in the algebra (no need of taking the real part) and reads In a similar fashion, it could be possible to show that also are legitimate expressions for the drift. All these new formulae are numerically different from those in Eqs. (29) and (33) but lead to the same results on average; clearly the advantage is that only half of the Lie derivative has to be computed.

C Fourier transforms with twisted boundary conditions
If f (x) is a periodic function defined on the L 4 lattice, its Fourier transform and inverse are where p is the quantised vector p = 2π L (n 1 , n 2 , n 3 , n 4 ) and the sum is to be read p = L−1 n 1 ,n 2 ,n 3 ,n 4 =0 . Antiperiodicity in the directionν would lead again to Eq. (72) but with a quantised momentum ( p ) μ = 2π L n μ + π L δ μν . Twisted boundary conditions on a plane Let us consider some N c × N c matrix M(x) (which for example can be a gauge link or the vector potential seen as matrices in colour space, or a fundamental fermion field seen as a matrix in colour-smell). We impose twisted boundary condition in thê 1,2 plane so that with 2 1 = z 1 2 , z = z 12 ∈ Z N . If we had just (anti)periodic boundary conditions, we would treat the matrix as N 2 c independent scalar functions; twisted boundary conditions actually couple the different components, therefore in order to expand M(x) in plane waves we need to find a good basis for the matrix space: it can be proved (see Refs. [32,33]) that the Fourier transform and its inverse are where p = p + p ⊥ , p ⊥ is the quantised vector p ⊥ = 2π N c L (ñ 1 ,ñ 2 , 0, 0) and the sum is to be read p ⊥ = N c −1 n 1 ,ñ 2 =0 . The matrices p ⊥ form the sought basis in the matrix space: assuming a twist with z = exp(2πi/N c ), we can choose for example A different choice for z would have somehow reshuffled the exponents in Eq. (75). We see that the Fourier transform of M(x) is a scalar functionM( p ) p ⊥ , but momentum has a finer resolution compared to (anti)periodic boundary conditions: spatial and colour degrees of freedom mix in momentum space. Moreover, traceless matrices naturally do not have a zero momentum component, becausẽ

Twisted boundary conditions in three directions
The conditions in Eq. (73) are supplemented by and then to each of these we apply the FFT, At the end, N 2 c projections and N 2 c FFTs have been performed. The inverse transform will be simplŷ followed by Note thatM( p ) p ⊥ is a scalar function but the dependence on p ⊥ is throughñ 1 ,ñ 2 , where each integer runs from 0 to N c − 1: this allows a representation of the Fourier transform again with a N c × N c matrix field, M( p ) ñ 1ñ2 . Of course this has to be understood only as a useful representation of the momentum degrees of freedom, not as a matrix in colour space.

D Autocorrelations and cross-correlations
We consider a sample {a i , b i } N i=1 of measures of two observables A, B taken from the stochastic process at equilibrium. Let A = a, B = b be the expectation values respectively of the observables A, B.The cross-correlation function is defined as where we used the fact that the expectation value is not dependent on i because the equilibrium distribution is timeindependent. The cross-correlation function is not an even function, AB (−t) = B A (t). In particular, AB (0) = Cov(A,B) is the covariance between A and B. The averagē a = 1 N N i=1 a i is a stochastic variable that satisfies ā = a. The covariance between the estimatorsā andb is but since the cross-correlation function is expected to drop exponentially at large times, it is possible to approximate with the integrated cross-correlation time We expect τ int AB = 1 2 when the observable B has some dependence on A. If B is independent of A, we can assume τ int AB = 1 2 . An estimator for the cross-correlation function is

E Twisted lattice perturbation theory
Twisted lattice perturbation theory for the a pure gauge theory was introduced in Ref. [32] (see also Ref. [66]). Recently, the computation of Wilson loops has been treated in detail in Ref. [48]. Here we focus on two vertices, introducing Wilson and staggered fermions with smell in the fundamental representation. Feynman rules are fairly similar to those of lattice perturbation theory, apart from phases in propagators and vertices; all phases cancel in the first-order computations we considered. We recall also that the sum over momenta is inherited from the Fourier transform in Appendix C, and each fermion loop has to be divided by N c , i.e. by the numbers of smells running in the loop. The function f ( p ⊥ , p ⊥ ) = z −ñ 1ñ 2 is introduced for convenience. The gluon propagator is where ξ is the gauge fixing parameter; note that the traceless property of the gauge field forces the propagator to vanish for p ⊥ = 0. The Wilson and staggered propagators are defined respectively in Eqs. (30) and (35). Below we write the fermion-fermion-gluon and fermion-fermiongluon-gluon vertices in the Wilson and staggered case; p 1 , p 2 are respectively the incoming and outgoing momenta of the fermions, k 1 , k 2 are the outgoing momenta of the gluons.
Wilson fermions

Staggered fermions
Here momentum conservation is made explicit, because the vertices are not diagonal in momentum space.
V f f g ( p 1 , p 2 , k 1 ) μ = −ig f (k 1⊥ , p 2⊥ ) cos p 2 + k 1 2 μ Fig. 15 Sketch of the PRlgt auxiliary links in a plane of a 2 × 2 lattice. Physical links are in black and sites identified by the same symbol represent the same physical site. Dashed lines highlight links that are allocated but do not participate in the update. In the left panel direction 1 is twisted,2 is not: in red there are the auxiliary sites (two forward, three backward) and the red links beginning there correspond to phys-ical links twisted according to the matrix 1 . In the right panel botĥ 1 and2 are twisted directions: in blue there are auxiliary sites whose links are twisted according to the matrix 2 . In the latter case there are sites which pass the boundary in two twisted directions: the green links undergo both the 1 and the 2 twist (the two operations commute by definition) × sin p 2 + k 1 2 + k 2 2 μ × δ μνδ (4) (− p 1 + k 1 + k 2 + p 2 + πμ)δ k 1⊥ − p 1⊥ +k 2⊥ + p 2⊥ ,0 (89b)

F Code development for NSPT
We developed two independent NSPT codes in order to crosscheck and improve our implementation. PRlgt 10 stems from the first NSPT codes developed by the Parma lattice gauge theory group, allowing for SU(3) simulations with Wilson fermions. We implemented twisted boundary conditions, smell for Wilson fermions and added support for SU (2) simulations. The code is tailored for perturbation theory. The underlying idea is to have base classes ptSU2 and ptSU3 that describe perturbative matrices. The operator * is overloaded with the Cauchy product, so that it is possible to write the product of two series in a natural way. This is one of the operations that, especially at high orders, becomes very time-consuming: thus, having perturbative matrices as base classes allows to keep the perturbative orders close in memory and to speed up the multiplication of series. In par- 10 For recent developments on the code see Ref. [67]. ticular, the perturbative expansion is hardcoded to start from 1 for an element of the group and from 0 for an element of the algebra, in order to avoid multiplying by the identity or zero matrix; this choice also improves numerical stability in keeping the series within the group or algebra. All the other structures are built from the base classes by adding Lorentz, Dirac or lattice degrees of freedom. The fermion field too is described by matrices in colour-smell space. The update of the configuration is done one link at a time: this is possible, faster and less memory consuming for the first order integrator we are using; indeed the staples around a link can be computed also if the neighbour links have already been updated, since the effect of doing so gives higher-order effects in the time step. Twisted boundary conditions are implemented ad hoc for the Wilson action, as shown in Fig. 15: a system of twisted copies of the links on the boundary is updated at each Langevin step. The code makes heavy use of multithreading in all loops over lattice sites. Even though the performance of PRlgt is extremely good for small lattices, it is hard to scale to large volumes due to the scalar nature of the code.
We have also developed the GridNSPT code, 11 based on the Grid library [68]. GridNSPT has been debugged against PRlgt, and we are able to obtain the very same outputs from these two completely different implementations (but staggered fermions have been implemented in GridNSPT only). The Grid library provides an environment where mes-sage passing, multithreading and vector parallelism are fully exploited: the lattice is geometrically decomposed into MPI domain, each one mapped to a set of processors; it is also overdecomposed over virtual nodes in order to fill a SIMD vector, assuring very high vectorisation efficiency. For example, on KNL and Skylake machines we can exploit the AVX-512 instruction set and a SIMD vector has room for 4 complex numbers in double precision; the virtual node decomposition results in the layout 1.1.2.2, where we are referring respectively to the coordinates x.y.z.t. Within the MPI task, multithreading is automatic because it is included in the closure of Grid lattice object expression templates. Grid incorporates C++11 internal template classes representing scalars, vector or matrices. We introduced a new template class representing a perturbative series, that embeds the overloading of the * operator.
template<class vtype, int Np> class iPert { vtype _internal[Np]; }; All the structures are tensors built from these templates: for example, the gauge field is Lattice<iVector<iScalar <iPert<iMatrix<vComplexD,Nc>,Np>>,Nd>>, where (starting from the outer template) we have the lattice, Lorentz, spin, perturbative, colour structure and the base type is a vectorised complex number in double precision. With this in place, every operation in Grid is performed consistently with almost no modification. We rely on the Grid library for the optimal implementation of the gauge action and for the Wilson and staggered fermion kernel. Twisted boundary conditions have been implemented modifying the covariant circular shifts. Even though GridNSPT lacks of many optimisations compared to PRlgt (for example the Langevin update is not performed one link at a time, but all operations and shifts are performed on the lattice as whole), it allows to have a more flexible environment and to scale easily end very efficiently to large volumes.

G The nearest covariance matrix
If C is a covariance matrix, the corresponding correlation matrix is defined aŝ where S is the matrix which is equal to C on the diagonal and vanishes everywhere else.Ĉ has 1 on the diagonal by construction; it might have some negative or zero eigenvalue if the estimator used in the determination of the covariance does not guarantee positive definiteness. Given C, Higham's algorithm [69] allows to find the nearest (in a weighted Frobenius norm) positive semidefinite matrix with unit diagonal. The core of the procedure is alternating a projection P S onto the space of positive semidefinite matrices and a projection P U onto the matrices with unit diagonal. The projection P S (X ) = Y consists in • diagonalising X = U T U , where U is an orthogonal matrix and is a diagonal matrix with the eigenvalues of X on the diagonal • setting to zero all the negative elements in , obtaining • returning Y = U T˜ U .
The projection P U (X ) consists simply in putting 1 on the diagonal of X . We refer to the original work for the presentation and proof of the complete algorithm: after some iterations, the algorithm converges and returns a matrixĈ H which is positive semidefinite and has 1 on the diagonal.
However, the algorithm allowsĈ H to have some zero (within machine precision) eigenvalue, preventing the inversion of the covariance matrix. If this is the case, we additionally projectĈ H onto the space of positive definite matrices. This projection consists in where V is an orthogonal matrix and is a diagonal matrix with the eigenvalues ofĈ H on the diagonal • identifying = δλ max , where λ max is the maximum eigenvalue and δ is the tolerance of the projection • setting to all the diagonal elements of whose absolute value is smaller than , obtaining˜ • returningĈ P = V T˜ V .
In conclusion, the nearest covariance matrix is