SMD-based numerical stochastic perturbation theory

The viability of a variant of numerical stochastic perturbation theory, where the Langevin equation is replaced by the SMD algorithm, is examined. In particular, the convergence of the process to a unique stationary state is rigorously established and the use of higher-order symplectic integration schemes is shown to be highly profitable in this context. For illustration, the gradient-flow coupling in finite volume with Schr\"odinger functional boundary conditions is computed to two-loop (i.e. NNL) order in the SU(3) gauge theory. The scaling behaviour of the algorithm turns out to be rather favourable in this case, which allows the computations to be driven close to the continuum limit.


Introduction
Numerical stochastic perturbation theory (NSPT) [1][2][3] is a powerful tool that allows many interesting calculations in QCD and other quantum field theories to be performed to high order in the interactions. For technical reasons, the computations proceed in the framework of lattice field theory, but results for renormalized quantities in the continuum theory can then be obtained through an extrapolation to vanishing lattice spacing. NSPT can be highly automated and the application of the method in finite volume and to correlation functions of complicated composite fields gives rise to hardly any additional difficulties. Reliable extrapolations to the continuum limit require accurate data at several lattice spacings in the scaling region. NSPT calculations can therefore rapidly become large-scale projects, where computational efficiency is all-important. Traditionally, NSPT is based on the Langevin equation, but the success of the HMC algorithm [4] in lattice QCD suggests that the inclusion of a molecular-dynamics update step in the underlying stochastic process might be beneficial. Smaller autocorrelation times and an improved scaling behaviour towards the continuum limit could perhaps be achieved in this way. Moreover, through the use of highly efficient symplectic integration schemes, the systematic errors deriving from the discretization of the simulation time may conceivably be reduced.
NSPT based on the SMD (stochastic molecular dynamics, or generalized HMC) algorithm [5] has recently been briefly looked at in ref. [6] and was found to perform well. Here we establish the convergence of the algorithm to a unique stationary state and study its efficiency in the case of the gradient-flow coupling in the SU(3) gauge theory. Various technical problems are addressed along the way, among them the modifications required to ensure that the stochastic process does not run away in the gauge directions.

Stochastic molecular dynamics
In order to bring out the basic structure of the SMD-variant of NSPT most clearly, a generic system described by a set q = (q 1 , . . . , q n ) of real coordinates and an action S(q) is considered in this and the following two sections.

Preliminaries
The action S(q) is assumed to be differentiable and to have an expansion in powers of a coupling g of the form where S r (q) is a polynomial in q of degree d r ≥ 2. Moreover, it is taken for granted that the leading-order term S 0 (q) = 1 2 (q, ∆q) = 1 2 n k,l=1 is a strictly positive quadratic form in the coordinates. The observables O(q) of interest are assumed to be similarly expandable and their expectation values then have a well-defined perturbation expansion with coefficients given by Feynman diagrams as usual.

SMD algorithm
The SMD algorithm operates in the phase space of the theory and thus updates both the coordinates q and their canonical momenta p = (p 1 , . . . , p n ). An SMD update cycle consists of a momentum rotation followed by a molecular-dynamics evolution and, optionally, an acceptance-rejection step. The momenta are rotated in a random direction according to p → c 1 p + c 2 υ, (2.4) where the momentum υ is randomly chosen from a Gaussian distribution with mean zero and unit variance. The coefficients depend on the simulation time step ǫ > 0 and a parameter γ > 0 that controls the rotation angle.
In the second step, the molecular-dynamics equations ∂ t p = −∇S(q), ∂ t q = p, (2.6) are integrated from the current simulation time t to t+ǫ using a reversible symplectic integration scheme (see subsect. 2.3). The algorithm (momentum rotation followed by the molecular-dynamics evolution) simulates the canonical distribution 1 Z H e −H(p,q) , H(p, q) = 1 2 (p, p) + S(q), (2.7) provided the integration errors are negligible or an acceptance-rejection step is included in the update cycle which corrects for these [7]. Stochastic estimates of the expectation values (2.3) of the observables of interest are then obtained by averaging their values over a range of simulation time.

Integration schemes
The molecular-dynamics equations may be integrated by applying a sequence of the elementary steps I p,h : p → p − h∇S(q), (2.8) I q,h : q → q + hp, (2.9) to the current momenta and coordinates, with step sizes h proportional to ǫ. A wellknown example is the leapfrog integrator I p,ǫ/2 I q,ǫ I p,ǫ/2 , and several highly efficient schemes are described in ref. [8].
Integrators I ǫ of this kind are symplectic and they can be (and are here) required to be reversible, i.e. to be such that where P stands for the momentum reflection p → −p.

Stochastic perturbation theory
Stochastic perturbation theory [9,10] is usually derived from the Langevin equation by expanding the stochastic variables and the driving forces in powers of the coupling. In this section, another (although probably closely related) form of stochastic perturbation theory is discussed, which is obtained by expanding the SMD algorithm in the same way.

SMD algorithm at weak coupling
Since the acceptance-rejection step is not smooth in the coupling, its effects would be difficult to take into account in perturbation theory. In the following, the acceptancerejection step is therefore omitted, without further notice, and one is thus left with an algorithm that simulates the system only up to integration errors. The histories p(t), q(t) of the momenta and coordinates generated by the SMD algorithm depend on the coupling g through the force term in the integration step (2.8). In particular, they are smooth functions of the coupling and may consequently be expanded in the asymptotic series where the leading-order historiesp 0 (t),q 0 (t) coincide with the ones generated by the algorithm in the free theory with action (2.2). In terms of the coefficientsp r ,q r , the momentum rotation (2.4) becomeŝ and the molecular-dynamics integration steps (2.8),(2.9) assume the form The forcesF r in eq. (3.3) are given by and it is understood that all momenta and all coordinates are updated alternately so that the variables on the right of eqs. (3.3), (3.4) are always the current ones. For any given initial data, these rules completely determine the historiesp r (t),q r (t).

Perturbation expansion of observables
Similarly to the gradient of the action, any observable may be expanded in powers of the coupling. The coefficients k r (O) in the perturbation expansion of its expectation value (2.3) then coincide with the averages ofÔ r (q 0 (t), . . . ,q r (t)) over the simulation time t up to statistical (and integration) errors.

Convergence to a stationary state
Stochastic processes can run away or do not converge to a stationary distribution for other reasons. In the case of the stochastic perturbation theory described in sect. 3, the asymptotic stationarity of the underlying process can be rigorously shown if the simulation step size ǫ is sufficiently small. The range of step sizes, where convergence is guaranteed, depends on the chosen integration scheme for the molecular-dynamics equations and the matrix ∆ in the leading-order part (2.2) of the action.

Molecular-dynamics evolution in the free theory
If the coupling g is turned off, the molecular-dynamics equations become linear and their (approximate) integration from time t to t + ǫ amounts to a linear transformation of the current momenta and coordinates. The 2n × 2n matrix M in this equation has a block structure, with n × n blocks that are polynomials in ǫ and the matrix ∆ with some numerical coefficients. In particular, they are commuting real symmetric matrices. In the case of the leapfrog integrator, and explicit expressions for the blocks can be obtained for other popular integrators as well (see appendix A). The symplecticity and reversibility of the chosen integration scheme imply It follows from these relations that the Hamilton function H(p, q) = 1 2 (p, p) + 1 2 (q,∆q), (4.6) is exactly conserved by the integrator. In the following, ǫ and ∆ are assumed to be such that M qp is non-singular and∆ positive definite. Both conditions are met in the case of the leapfrog integrator if ǫ 2 ∆ < 4, where ∆ is the largest eigenvalue of ∆. The probability distribution is then well-behaved and preserved by the SMD algorithm, since it is preserved both by the momentum rotation and the molecular-dynamics evolution.

Convergence of the leading-order process
For any initial distribution of the momenta and coordinates, the SMD algorithm produces a sequence of distributions, which converges to the stationary distribution (4.8) in the free theory. One can show this by working out the action of an SMD update cycle on a given distribution, but the convergence of the algorithm may be established more easily starting from the identity where υ(u) is the momentum chosen randomly in the momentum rotation (2.4) at simulation time u = 0, ǫ, 2ǫ, . . . and The long-time behaviour of the momenta and coordinates thus depends on the properties of the matrixM . The blocks M pp , M pq , M qp and∆ are commuting real symmetric matrices. Since and since∆ is assumed to be positive definite, the eigenvalues of M pp have magnitude strictly less than 1. The eigenvalues ofM are then where µ runs through the eigenvalues of M pp . In particular, |λ ± | < 1 andM is thus a contraction matrix. The first term on the right of eq. (4.9) consequently dies away exponentially with increasing simulation time t. Since the random momenta are normally distributed, the momenta and coordinates at large t are then normally distributed as well, with mean zero and variances equal to their two-point autocorrelation functions. In the large time limit, the qq autocorrelation function, for example, is given by 14) and the other two-point functions by the pp, pq and qp blocks of the matrix on the right of eq. (4.13). The kernel K satisfies i.e. an inhomogeneous linear equation that has a unique solution, sinceM is a contraction matrix. A few lines of algebra then show that solves the equation. In particular, the variances of the momenta and coordinates at any fixed simulation time coincide with the ones of the stationary distribution (4.8), which proves that the latter coincides with the distribution simulated by the SMD algorithm.

Convergence beyond the leading order
The assumed structure of the action S(q) implies that the force in the moleculardynamics integration step (3.3) is of the form F r (q 0 , . . . ,q r ) = ∆q r +F ′ r (q 0 , . . . ,q r−1 ), (4.17) where the second term on the right is a polynomial in the coordinates up to order r − 1. If the history of the latter and the associated momenta is already known (including their values at the intermediate integration steps), the integration of the molecular-dynamics equations at order r thus amounts to solving an inhomogeneous linear recursion. A moment of thought then reveals that for all r ≥ 1. The "vertices"V r (p 0 ,q 0 ; . . . ;p r−1 ,q r−1 ) in this formula are polynomials in their arguments, whose exact form depends on both the integration scheme and the forces (4.17). Recalling eq. (4.9) and the fact thatM is a contraction matrix, the convergence of the autocorrelation functions of the momenta and coordinates at large times t may now be shown recursively from order 0 to any finite order r. Equation (4.18) actually allows the highest-order variables in any correlation function to be expressed through lower-order ones up to an exponentially decaying contribution. Clearly, in the largetime limit, the autocorrelation functions do not depend on the initial distribution of the variables and are stationary, i.e. invariant under time translations.
The expectation values of the coefficients in the perturbation expansion (3.6) of the observables coincide with a sum of autocorrelation functions of the coordinates at equal times. Their convergence at large times is therefore guaranteed as well.

Summary
The discussion in this section shows that the SMD algorithm converges to all orders of perturbation theory if the matrix ∆ is strictly positive and if ǫ 2 ∆ < κ, where κ depends on the molecular-dynamics integrator. In the case of the leapfrog, the 2nd-order OMF and the 4th-order OMF integrators, κ is equal to 4, 6.51 and 9.87, respectively (cf. appendix A).

Stochastic perturbation theory in lattice QCD
With respect to the generic system considered so far, the situation in lattice QCD is complicated by the gauge symmetry and the quark fields. In this section, stochastic perturbation theory is first set up for the pure SU(N ) gauge theory. The modifications required for the damping of the gauge modes are then discussed and the section ends with a brief description of how the quarks can be included in the simulations.

Lattice fields
The lattice theory is set up on a T × L 3 lattice with periodic boundary conditions in the space directions and Schrödinger functional (SF) [11,12] or open-SF [13] boundary conditions in the time direction. In both cases the link variables U (x, µ) satisfy at time x 0 = T and additionally at x 0 = 0 if SF boundary conditions are chosen. The notation used in the following coincides with the one previously employed in ref. [13]. In particular, all dimensionful quantities are expressed in units of the lattice spacing.
The momenta π(x, µ) of the link variables U (x, µ) take values in the Lie algebra su(N ) of the gauge group. They vanish at the boundaries of the lattice, where the link variables are frozen to unity. The scalar product of any two momentum fields, includes a conventional weight factor which will reappear below in various expressions. Infinitesimal gauge transformations are fields ω(x) of su(N ) elements defined on the sites x of the lattice. They vanish at x 0 = T and must, furthermore, be constant at x 0 = 0 if SF boundary conditions are imposed. In this latter case, the fields may be split in two parts according to whereω(x) vanishes at x 0 = 0, T and is otherwise unconstrained. A possible choice of scalar product is then with w x = 2 at x 0 = 0, T and w x = 1 elsewhere. Independently of the boundary conditions imposed on the gauge field, the quark fields are required to vanish at both x 0 = 0 and x 0 = T . No weight factor is included in the scalar product of these fields.

Basic stochastic process
In the absence of quark fields, the SMD algorithm proceeds along the lines of sect. 2. For the action S(U ) of the gauge field any of the frequently used ones may be taken. The Hamilton function from which the SMD algorithm derives is then given by In particular, the random field in the momentum rotation must be distributed with probability density proportional to e − 1 2 (υ,υ) . If the weight factor w x,µ in the scalar product (5.2) is also included in the symplectic structure (i.e. the Poisson bracket), the molecular-dynamics equations assume the form For later convenience, the expressions on the right of the equations have been scaled by the bare gauge coupling g 0 , an operation that could be undone by rescaling the simulation time. The integration steps (2.8) and (2.9) are then given by Since it will be omitted in perturbation theory, the acceptance-rejection step needed to correct for the integration errors is not described here.

Perturbation expansion
In perturbation theory, the link variables are represented by a gauge potential A µ (x) through The gauge potential takes values in su(N ) and satisfies the same boundary conditions as the momentum field When the gauge and momentum fields are replaced by these expansions, the SMD algorithm leads to a hierarchy of stochastic equations as in the case of the generic system considered in sect. 3. At lowest order in the coupling, the momentum is rotated according to eq. (5.8) (withπ 0 in place of π) and the integration steps (5.10),(5.11) become where ∆ is the symmetric linear operator in the leading-order expression for the gauge action.

Damping of the gauge modes
The space H 1 of gauge potentials may be split into the subspace H L 1 of gauge modes and its orthogonal complement H T 1 . There is a one-to-one correspondence between the infinitesimal gauge transformations, introduced in subsect. 5.1, and the gauge modes through the forward difference operator which maps any field ν in the space H 0 of infinitesimal gauge transformations to a field dν ∈ H L 1 . The adjoint operator d * goes in the opposite direction and is defined by the requirement that for all ν ∈ H 0 and A ∈ H 1 . In particular, the subspace H T 1 coincides with the space of gauge potentials A satisfying d * A = 0 (d * is given explicitly in appendix B).
Since the gauge modes are annihilated by ∆, the H L 1 -component of the leadingorder gauge fieldÛ 0 performs a random walk and thus slowly runs away in the course of a simulation. Stability can be regained by applying a gauge transformation 20) to the fields after each SMD update cycle, where ω ∈ H 0 is a new field that is evolved together with the other fields. There are two update steps for this field, the first, being applied together with the momentum rotation and the second, at the end of the molecular-dynamics evolution. The parameter λ 0 > 0 controls the feedback from the current gauge field to the gauge-damping field ω(x) and can, in principle, be set to any value. Clearly, the history of gauge-invariant observables is not affected by these modifications of the SMD algorithm, but as discussed below (in subsect. 5.5), they have the desired damping effect on the gauge modes †.
Like the other fields, the gauge-damping field is expanded in a series in perturbation theory. At leading order, the new update steps are then 26) † In the continuous-time limit ǫ → 0, the time-dependence of the gauge-damping field is governed by a first-order differential equation and the modified algorithm integrates a form of the stochastic molecular-dynamics equations, which coincides with the standard one up to a time-dependent gauge transformation (see appendix C). 27) and the higher-order rules have the familiar hierarchical structure.

Long-time stationarity of the process
Although the stochastic process now includes variables without associated momenta, the discussion in sect. 4 applies here too with only minor changes. In particular, the convergence of the process to a unique stationary state is guaranteed to all orders of the coupling, if the matrixM describing the evolution of the fields (π 0 ,Û 0 ,ω 0 ) at leading order is a contraction matrix.
Since the gauge modes are zero modes of ∆, the subspace of field vectors is invariant under the action ofM . When restricted to its orthogonal complement (which is invariant too), the matrix describes the evolution of the H T 1 -components of π 0 ,Û 0 . The chosen boundary conditions (SF or open-SF) imply the strict positivity of ∆ in this subspace andM is therefore a contraction matrix there if where κ depends on the molecular-dynamics integrator (cf. subsect. 4.4).
In the subspace (5.28) of the gauge modes, on the other hand, the action ofM amounts to applying the linear transformation to the fields. The associated eigenvalues are equal to c 1 or to where µ 2 is any eigenvalue of d * d. This operator has no zero modes and some simple estimates then show thatM is a contraction matrix in the subspace of gauge modes, provided the bound holds.
Convergence of the stochastic process to equilibrium is thus guaranteed if the inequalities (5.29) and (5.32) are both satisfied. With the chosen boundary conditions, where k is equal to 1, 5/3 and 3.648 in the case of the Wilson [14], tree-level Symanzik improved [15,16] and Iwasaki [17] gauge action, respectively. In practice, the constraints (5.29) and (5.32) therefore tend to be fairly mild and the main concern is to ensure that the integration errors are sufficiently small at the chosen values of ǫ.

Inclusion of the quark fields
As in the case of the HMC algorithm [4], the effects of the quarks can be included in SMD simulations by adding a multiplet of pseudo-fermion fields to the theory with the appropriate action. The details are not important here and it suffices to know that the action is a sum of terms, one per pseudo-fermion field φ, of the form with R(U ) being some gauge-covariant linear operator (see ref. [18], for example).
Only few modifications of the SMD update cycle are required in the presence of the pseudo-fermion fields. Similarly to the momentum field, each field is rotated according to † at the beginning of the cycle, where η is a Gaussian random field with mean zero and unit variance. The molecular-dynamics evolution of the gauge and momentum field then proceeds at fixed pseudo-fermion fields, with the contribution of the pseudofermion action to the driving force properly taken into account. Clearly, the gauge transformation (5.19), (5.20) applied at the end of the update cycle must be extended to the pseudo-fermion fields. When the algorithm is expanded in powers of the coupling g 0 , the renormalization of the quark masses should be taken into account so that the masses in the leadingorder stochastic equations are the renormalized ones. The pseudo-fermion fields decouple from the other fields at lowest order and are simulated exactly by the random rotation (5.35). Convergence of the stochastic process to a unique stationary state is then again guaranteed to all orders, if the bounds (5.29) and (5.32) are satisfied. † There is no reason other than simplicity to set the parameter γ that determines the rotation angle to the same value for all fields.

Computation of the gradient-flow coupling
The gradient-flow coupling in finite volume with SF boundary conditions has recently been used in step scaling studies of three-flavour QCD [20,21]. Such calculations serve to relate the low-energy properties of the theory to the high-energy regime, where contact with the standard QCD parameters and matrix elements can be made in perturbation theory [22] (see ref. [23] for a review).
In the following, the perturbation expansion of the gradient-flow coupling is determined in the SU(3) Yang-Mills theory up to two-loop order, using the SMD-variant of NSPT described in the previous section. To one-loop order, the expansion coefficient in the MS scheme of dimensional regularization is known in infinite volume since a while [24], but a huge effort plus the best currently available techniques were required to be able to extend this calculation to the next order [25]. In finite volume with SF boundary conditions, these techniques do not apply and a similar analytical computation may be practically infeasible at present.

Definition of the coupling
The renormalized coupling considered in this paper belongs to a family of couplings based on the Yang-Mills gradient flow. Explicitly, it is given by [19] where E(t, x) denotes the Yang-Mills action density at flow time t and position x, c is a parameter of the scheme and the proportionality constant k is determined by the requirement thatḡ 2 coincides with g 2 0 to lowest order of perturbation theory. Most of the time c will be set to 0.3, which implies a localization range of the action density of about 0.3 × L.
On the lattice, the gradient flow is implemented as in ref. [13]. The Wilson action, with boundary terms so as to ensure the absence of O(a) lattice effects in the flow equation, is thus used to generate the flow. For the action density E(t, x) in eq. (6.1) either the Wilson action density or the square of the familiar symmetric "clover" expression for the gauge-field tensor is inserted. Furthermore, alternative couplings, where the full action density is replaced by its spatial or time-like part, are considered as well.
All in all this makes 6 different action densities and couplings, labeled w, ws, wt, c, cs, and ct, where the letters stand for Wilson, clover, space and time, respectively (E ct , for example, denotes the clover expression for the time-like part of the action density).

Expansion in powers of α s
The gradient flow coupling has an expansion in powers of the running coupling α s in the MS scheme of dimensional regularization at momentum scale q, with coefficients k 1 , k 2 , . . . depending on q and L. If q is scaled with L like the coefficients are of order 1 and independent of L in the continuum theory [26]. In the following k 1 and k 2 are computed by combining the expansion obtained in NSPT with the relation between α s and the bare coupling. For the Wilson gauge action (which is the action used in NSPT), the coefficients d 1 and d 2 are accurately known [27]. The coefficients k 1 , k 2 determined in this way depend on the spacing of the simulated lattices so that an extrapolation to the continuum limit is then still required.

Computation of the coefficients E k in NSPT
In order to cancel the O(a) lattice effects in E k , the action of the theory must include boundary counterterms at x 0 = 0 and x 0 = T with a tuned coupling [11] † as was shown long ago [28][29][30]. The counterterms lead to further terms in the forces that drive the molecular-dynamics evolution, but do not require a modification of the † In ref. [13] the improvement coefficient c t is denoted by c ′ G . general framework described in sect. 5. Alternatively, the effects of the counterterms can be included in the calculations by treating c t − 1 as a second coupling. Once a representative sample of gauge configurations (5.12) has been generated, the stochastic estimation of the gradient-flow coupling requires E(t, x) to be expanded in powers of g 0 for each of these configurations. To this end, the gauge field V t (x, µ) at flow time t is represented by a gauge potential B µ (t, x) according to The numerical integration of the flow equation, using a Runge-Kutta integrator, for example, then amounts to applying a sequence of integration steps to the expansion coefficients of the field as in the case of the integration of the molecular-dynamics equations. Gauge damping is however not required here, since the linearized flow is transversal and leaves the gauge modes unchanged.

Simulation parameters and tables of results
The parameters of the NSPT runs performed for the "measurement" of E 0 , E 1 , E 2 are listed in table 1. In all cases, the gauge-damping parameter λ 0 was set to 2 and the 4th-order OMF integrator was employed for the molecular-dynamics equations (cf. appendix A). Measurements were made after every ∆t ms /ǫ update cycles using a 3rd-order Runge-Kutta integrator for the gradient flow [24], with step size varying from 0.002 at small flow times to 0.1 at large times. With this scheme, the gradientflow integration errors are guaranteed to be completely negligible with respect to the statistical errors. The number N ms of measurements made is listed in the last column of table 1 (the programs that were used in these simulations can be downloaded from http://luscher.web.cern.ch/luscher/NSPT).
At the chosen values of the parameters, the bounds (5.29) and (5.32) are satisfied by a wide margin so that the convergence of the SMD algorithm is rigorously guaranteed. The results obtained on the simulated lattices for the expansion coefficients k 1 and k 2 and their statistical errors are listed in appendix D, for c ∈ {0.2, 0.3, 0.4} and all choices w,c,ws,. . . of the action density (cf. subsect. 6.1).

Statistical and systematic errors
The values of k 1 and k 2 obtained in NSPT depend on the scheme parameter c, the lattice size L (in lattice units), the simulation step size ǫ and the SMD parameter γ. No attempt is made here to determine all these dependencies in detail. Instead some basic facts and empirical results are discussed that helped controlling the errors in the case of the simulations reported in this paper (see refs. [6,31] for related complementary studies of the φ 4 theory).

Autocorrelations and statistical variances
Usually the variances of the observables are a property of the simulated field theory and hence independent of the simulation algorithm. In NSPT this is not so, because the algorithm is not exact, but mainly because the square of the order-r coefficient O r of an observable O in general does not coincide with the order-2r coefficient of another observable. The time average ofÔ 2 r and thus the variance ofÔ r are then not necessarily determined by the theory alone. For illustration, the dependence on γ of the variances of the coefficients in (where x 0 = L/2 and √ 8t = cL as before) is shown in fig. 1. As will be discussed in subsect. 7.3, the integration errors are negligible with respect to the rapid growth of the variances of the one-and two-loop coefficients seen at small γ, which is therefore entirely an effect of the change of algorithm.
In practice one would like to minimize the computational work required to obtain the calculated coefficients to a given statistical precision. The simulation algorithm should thus be tuned so as to minimize the products of the integrated autocorrelation times and variances of the order-r coefficientsÔ r of the observables O of interest. Empirical studies show that the two factors often work against each other, i.e. algorithms tuned for small autocorrelations tend to give large variances and vice versa. The autocorrelation times of the coefficientsÊ c k , for example, grow monotonically with γ (see table 2). At the chosen point in parameter space, the associated products (7.2) are then minimized at values of γ around 0.5, 1.0 and 2.0 for k = 0, 1 and 2. The example shows that there may be no uniformly best choice of γ, but a range of values, where all coefficients of interest are obtained reasonably efficiently.

Critical slowing down
The behaviour of the autocorrelation times and variances near the continuum limit L → ∞ depends on the simulation algorithm and the observables considered. When NSPT is based on the Langevin equation, the autocorrelations of the coefficients of multiplicatively renormalizable quantities can be shown to diverge proportionally to L 2 [32,33], while their variances grow at most logarithmically [34].
At large γ, the variant of NSPT studied here is expected to behave similarly, since the SMD algorithm then effectively integrates the Langevin equation. Choosing γ to Table 2. Autocorrelation times ofÊ c k and associated products ( depend on L in some particular way may, however, conceivably lead to an improved scaling behaviour. In the free theory, for example, the autocorrelation times grow proportionally to L rather than L 2 if γ is scaled like 1/L [5]. Beyond the leading order, the situation is complicated by the algorithm-dependence of the variances and the effects of the parameter tuning are then not easy to predict. Considering the fact that the computational cost of the measurements ofÊ k tends to be larger than the one of the SMD update cycles, the parameters of the runs on the large lattices (those with L = 24, . . . , 40 in table 1) were chosen so that subsequent measurements are practically uncorrelated. At fixed γ = 3 the required increase with L of the measurement time separation ∆t ms then turned out to be quite moderate. Moreover, from L = 24 to L = 40, the variances ofÊ k grow only slowly (at c = 0.3 and for k = 0, 1 and 2 by about 2, 19 and 30 percent).

Integration errors
As already noted in appendix A, the theory is very accurately simulated at leading order if the 4th-order OMF integrator is used for the molecular-dynamics equations.
The expectation values of the coefficientÊ 0 calculated in the runs listed in table 1 in fact all agree with the known analytic formula [19,13] within a relative statistical precision of about 2 × 10 −4 .
Beyond the leading order, the integration errors remain difficult to detect in empirical studies (see fig. 2). The stability bounds (5.29),(5.32) are respected in all these runs and the coefficientsĤ 0 , . . . ,Ĥ 4 of the Hamilton function H are accurately conserved, with deficits decreasing like ǫ 5 (as has to be the case in the asymptotic regime of a 4th-order integrator). It thus seems safe to conclude that the integration errors in the tests reported in fig. 2 are smaller than the statistical errors. Apart from adjustments to match the target statistics, the step sizes in the runs listed in table 1 were, for L ≥ 24, scaled proportionally to 1/L so that the integration errors fall off like 1/L 4 at large L and thus much more rapidly than the leading lattice effects. Since the statistical errors are approximately the same in all runs, the scaling of the step sizes may be a luxury, but was applied here as a safeguard measure against an enhancement of systematic errors through the continuum extrapolation.

Extrapolation to the continuum limit
With O(a)-improvement in place, the L-dependence of the one-loop coefficient k w 1 is asymptotically given by [35,36] where the leading-order term, a 0 , coincides with the coefficient k 1 in the continuum theory. The available data are well described by this asymptotic expression down to L = 10 if c ≥ 0.3 and L = 12 if c = 0.2 (see fig. 3). Note the fact that the data points in fig. 3 are uncorrelated and that random fluctuations by more than one standard deviation are consequently not improbable in a sample of 10 points. Plots of the other one-loop coefficients k c 1 , k ws 1 , . . . look much the same. Even though only a short extrapolation is required to reach the continuum limit, there is no way the extrapolation errors can be rigorously assessed. Following stan- dard practice, the errors are estimated by varying the fit procedure within reasonable bounds. In particular, fits are discarded if the quality of the fit is low or if the fit parameters are poorly determined. Another indication of the size of the extrapolation errors is provided by the deviations of the results obtained using the "w" and the "c" data, respectively. The errors of the results quoted in table 3 include both the statistical and an estimate of the extrapolation errors.
With increasing loop order, the continuum extrapolations tend to become more difficult, partly because one loses statistical precision and partly because the asymptotic formulae describing the leading lattice effects have more and more terms. In particular, with respect to the one-loop coefficients, the coefficients at the next order are obtained about 10 times less precisely and their asymptotic form includes an additional logarithm. Combinations of the logarithms in this expression only too easily accurately approximate a constant in the range of the available data, but may blow up closer to the continuum limit and then strongly influence the result of the extrapolation.
In the present context, the goal is not to determine the values of the coefficients a 1 , b 1 and c 1 , but to perform a short extrapolation of the data to the continuum limit. A sensible way to stabilize fits of the data by the asymptotic expression (7.4) is then to constrain the minimization of the likelihood function to the directions in parameter space orthogonal to its flattest directions (see fig. 4). The results quoted in table 3 were obtained in this way and by varying the fit procedure, as in the case of the one-loop coefficients, in order to assess the extrapolation errors.

Miscellaneous remarks
Lattice effects. Since the smoothing range of the gradient flow decreases with c, it is no surprise that the coefficients k 1 , k 2 calculated in NSPT are found to be increasingly sensitive to lattice effects when c is lowered. The continuum extrapolation of the data then becomes more and more difficult and eventually requires larger lattices to be simulated. Infinite-volume limit. The gradient-flow coupling in infinite volume runs with the flow time t and may be expanded in powers of α s at scale q = (8t) −1/2 , as in eq. (6.2), the one-and two-loop coefficients in the continuum limit being k 1 = 1.0978(1) [24] and k 2 = −0.9822(1) [25]. In the finite-volume scheme considered in this paper, and after passing to the continuum limit, the infinite-volume limit is reached by sending c to zero. The results listed in table 3 cannot be reliably extrapolated to c = 0, but the values of k 1 , . . . , k t 2 at c = 0.2 are actually already quite close to the infinite-volume values quoted above.
Three-loop β-function. The L-dependence of the gradient-flow coupling α =ḡ 2 /4π is governed by the renormalization group equation L ∂α ∂L = β 0 α 2 + β 1 α 3 + β 2 α 4 + . . . , (7.5) where β 0 = 11(2π) −1 and β 1 = 51(2π) −2 are the universal one-and two-loop coefficients of the Callan-Symanzik β-function. Using the results quoted in table 3, the three-loop coefficient may be calculated in a few lines (see table 4). The coefficient thus has opposite sign and is significantly larger than in the MS scheme, particularly so at c = 0.4 and if the spatial part of the Yang-Mills action density is inserted in the definition (6.1) of the coupling.

Conclusions
In stochastic perturbation theory the fields are represented by a series of coefficient fields that solve the underlying stochastic equation order by order in the interactions. Beyond the leading order, the probability distributions of the coefficient fields are however not a priori known and actually depend on the chosen stochastic process. The variances of the observables of interest must consequently be expected to vary  (8) with the parameters of the simulation algorithm, an effect that tends to considerably complicate the situation with respect to the one in simulations based on importance sampling.
The SMD algorithm provides a technically attractive basis for NSPT. Compared to the traditional version of NSPT [1][2][3], where the starting point is the Langevin equation, a significantly improved efficiency is achieved, particularly so near the continuum limit. Moreover, the available highly accurate integrators for the moleculardynamics equations allow the integration errors to be easily kept under control. For the reasons mentioned above, some tuning of the friction parameter γ is however required and must take into account the variances of the observables of interest.
The inclusion of the quark fields in the SMD algorithm along the lines of sect. 5 is straightforward and is not expected to slow down the simulations by a large factor [3]. In general, the cost of NSPT computations very much depends on the observables considered, the order in perturbation theory and the desired precision of the results.
The statistical errors of the expansion coefficients k l of the gradient-flow coupling, for example, appear to grow rapidly with the loop order l. In practice some loss of precision from one order to the next is tolerable, since the coefficients get multiplied by α l+1 s in the perturbation series (6.2). An extension of the computations to threeloop order would however only make sense if k 1 and k 2 are recalculated with errors about an order of magnitude smaller than the ones quoted in table 3. Furthermore, the relation between α s and the bare coupling would need to be worked out to three loops.
M.D.B. would like to thank Marco Garofalo, Dirk Hesse, and Tony Kennedy for a pleasant collaboration on related investigations, Alberto Ramos, Stefan Sint, and Rainer Sommer for valuable discussions and the Theoretical Physics Department at CERN for the kind hospitality extended to him. The computations reported in this paper were performed on the SuperMUC machine at the Leibniz Supercomputing Centre in Munich (project ID pr92ci) and dedicated HPC clusters at DESY-Zeuthen and at CERN. We thank all these institutions for the generous support given to this project. and the matrices M pq and M qp are polynomials in ǫ 2 ∆ of degree 5 and 4 times ǫ∆ and ǫ, respectively. In this case, M qp is non-singular and ∆ = ∆ 1 + a 1 ǫ 4 ∆ 2 + a 2 ǫ 6 ∆ 3 + . . . , (A.6) is positive if ǫ 2 ∆ < 9.87. As is evident from eqs. (A.6),(A.7), this integrator is impressively accurate. If ∆ = 16 and ǫ = 0.2, for example, the relative deviation of∆ from ∆ is at most 1.6 × 10 −5 and thus about a factor 12 smaller than the deviation in the case of the 2nd-order OMF integrator (with an adjusted step size of ǫ = 0.08).

Appendix B. Explicit form of d *
The operator d * depends on the chosen boundary conditions and the scalar products in the field spaces H 0 and H 1 (cf. subsect. 5.1). For SF boundary conditions while in the case of open-SF boundary conditions,

Appendix C. Gauge-damped stochastic equations
In the continuous-time limit ǫ → 0, the gauge-damped SMD algorithm described in sect. 5 integrates the stochastic equations where π t , U t , ω t denote the momentum, gauge and gauge-damping fields at simulation time t (and the link field C t is given by eq. (5.23) with U replaced by U t ). The normalization of the Gaussian random field is such that η a t (x, µ)η b s (y, ν) = 2γw x,µ δ ab δ µν δ(t − s)δ xy (C.4) and ∇ µ ω t (x) = U t (x, µ)ω t (x +μ)U t (x, µ) −1 − ω t (x) (C.5) stands for the gauge-covariant gradient of the gauge-damping field. As in the case of the SMD algorithm, the gauge-damping terms can be removed from the continuous-time process by applying a time-dependent gauge transformation to the fields π t , U t and η t . Up to a rescaling of the simulation time and the momentum field by the coupling g 0 , eqs. (C.1),(C.2) then reduce to the standard stochastic molecular-dynamics equations [5]. The gauge damping can also easily be shown to coincide with the standard one in the Langevin limit γ → ∞ [37].