SMD-based numerical stochastic perturbation theory

TL;DR

This paper explores a variant of numerical stochastic perturbation theory using the SMD algorithm, demonstrating its convergence, efficiency with higher-order integrators, and successful two-loop calculations in SU(3) gauge theory.

Contribution

It introduces and rigorously analyzes a new SMD-based approach to numerical stochastic perturbation theory, showing improved convergence and computational advantages.

Findings

Proves convergence to a unique stationary state.

Shows higher-order symplectic schemes are highly effective.

Performs two-loop calculations in SU(3) gauge theory successfully.

Abstract

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.

Equations

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad S(q)=\sum_{r=0}^{\infty}g^{r}S_{r}(q),\hfill\hbox to0.0pt{\hss{\tenrm(2.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad S_{0}(q)=\hbox{$1\over 2$}(q,\Delta q)=\hbox{$1\over 2$}\sum_{k,l=1}^{n}q_{k}\Delta_{kl}q_{l}\hfill\hbox to0.0pt{\hss{\tenrm(2.2)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\langle{\cal O}\rangle={1\over{\cal Z}_{S}}\int{\tenrm d}q_{1}\ldots{\tenrm d}q_{n}\,{\cal O}(q)\kern 0.5pt{\tenrm e}^{-S(q)}\hfill\hbox to0.0pt{\hss{\tenrm(2.3)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad p\to c_{1}p+c_{2}\upsilon,\hfill\hbox to0.0pt{\hss{\tenrm(2.4)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad c_{1}={\tenrm e}^{-\gamma\epsilon},\qquad c_{2}=(1-c_{1}^{2})^{1/2},\hfill\hbox to0.0pt{\hss{\tenrm(2.5)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\partial_{t}p=-\nabla S(q),\qquad\partial_{t}q=p,\hfill\hbox to0.0pt{\hss{\tenrm(2.6)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad{1\over{\cal Z}_{H}}{\tenrm e}^{-H(p,q)},\qquad H(p,q)=\hbox{$1\over 2$}(p,p)+S(q),\hfill\hbox to0.0pt{\hss{\tenrm(2.7)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad I_{p,h}:\;p\to p-h\nabla S(q),\hfill\hbox to0.0pt{\hss{\tenrm(2.8)}}\cr\vskip 8.61108pt\cr 0.0pt{\hfil$\displaystyle\qquad I_{q,h}:\;q\to q+hp,\hfill\hbox to0.0pt{\hss{\tenrm(2.9)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad I_{\epsilon}PI_{\epsilon}=P,\hfill\hbox to0.0pt{\hss{\tenrm(2.10)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad p(t)=\sum_{r=0}^{\infty}g^{r}\hat{p}_{r}(t),\qquad q(t)=\sum_{r=0}^{\infty}g^{r}\hat{q}_{r}(t),\hfill\hbox to0.0pt{\hss{\tenrm(3.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\hat{p}_{r}\to c_{1}\hat{p}_{r}+\delta_{r0}c_{2}\upsilon\hfill\hbox to0.0pt{\hss{\tenrm(3.2)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad I_{\hat{p}_{r},h}:\;\hat{p}_{r}\to\hat{p}_{r}-h\hat{F}_{r}(\hat{q}_{0},\ldots,\hat{q}_{r}),\hfill\hbox to0.0pt{\hss{\tenrm(3.3)}}\cr\vskip 10.76385pt\cr 0.0pt{\hfil$\displaystyle\qquad I_{\hat{q}_{r},h}:\;\hat{q}_{r}\to\hat{q}_{r}+h\hat{p}_{r}.\hfill\hbox to0.0pt{\hss{\tenrm(3.4)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\nabla S(q)=\sum_{r=0}^{\infty}g^{r}\hat{F}_{r}(\hat{q}_{0},\ldots,\hat{q}_{r})\hfill\hbox to0.0pt{\hss{\tenrm(3.5)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad{\cal O}(q)=\sum_{r=0}^{\infty}g^{r}\hat{{\cal O}}_{r}(\hat{q}_{0},\ldots,\hat{q}_{r})\hfill\hbox to0.0pt{\hss{\tenrm(3.6)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\langle{\cal O}\rangle=\sum_{r=0}^{\infty}k_{r}({\cal O})g^{r}\hfill\hbox to0.0pt{\hss{\tenrm(3.7)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\pmatrix{p\cr q\cr}\to M\pmatrix{p\cr q\cr}\hfill\hbox to0.0pt{\hss{\tenrm(4.1)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad M=\pmatrix{M_{pp}&M_{pq}\cr M_{qp}&M_{qq}\cr},\hfill\hbox to0.0pt{\hss{\tenrm(4.2)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad M_{pp}=M_{qq}=1-\hbox{$1\over 2$}\epsilon^{2}\Delta,\quad M_{pq}=-\epsilon\Delta\bigl{(}1-\hbox{$1\over 4$}\epsilon^{2}\Delta\bigr{)},\quad M_{qp}=\epsilon,\hfill\hbox to0.0pt{\hss{\tenrm(4.3)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad M_{pp}=M_{qq},\hfill\hbox to0.0pt{\hss{\tenrm(4.4)}}\cr\vskip 8.61108pt\cr 0.0pt{\hfil$\displaystyle\qquad M_{pp}M_{qq}-M_{pq}M_{qp}=1.\hfill\hbox to0.0pt{\hss{\tenrm(4.5)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\hat{H}(p,q)=\hbox{$1\over 2$}(p,p)+\hbox{$1\over 2$}(q,\hat{\Delta}q),\hfill\hbox to0.0pt{\hss{\tenrm(4.6)}}\cr\vskip 10.76385pt\cr 0.0pt{\hfil$\displaystyle\qquad\hat{\Delta}=-M_{pq}(M_{qp})^{-1}=\Delta+{\tenrm O}(\epsilon),\hfill\hbox to0.0pt{\hss{\tenrm(4.7)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\hat{P}(p,q)\propto{\tenrm e}^{-\hat{H}(p,q)}\hfill\hbox to0.0pt{\hss{\tenrm(4.8)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\pmatrix{p(t)\cr q(t)}=\tilde{M}^{t/\epsilon}\pmatrix{p(0)\cr q(0)}+c_{2}\sum_{u=0}^{t-\epsilon}\tilde{M}^{(t-u)/\epsilon-1}M\pmatrix{\upsilon(u)\cr 0},\hfill\hbox to0.0pt{\hss{\tenrm(4.9)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\tilde{M}=M\pmatrix{c_{1}&0\cr 0&1\cr}.\hfill\hbox to0.0pt{\hss{\tenrm(4.10)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad M_{pp}^{2}=1-M_{qp}\hat{\Delta}M_{qp}\hfill\hbox to0.0pt{\hss{\tenrm(4.11)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\lambda_{\pm}=\hbox{$1\over 2$}\bigl{\{}(1+c_{1})\mu\pm\sqrt{(1+c_{1})^{2}\mu^{2}-4c_{1}}\bigr{\}},\hfill\hbox to0.0pt{\hss{\tenrm(4.12)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\langle q(t)q(s)\rangle_{\upsilon}\mathrel{\mathop{=}_{t\geq s}}\bigl{\{}\tilde{M}^{(t-s)/\epsilon}K\bigr{\}}_{qq},\hfill\hbox to0.0pt{\hss{\tenrm(4.13)}}\cr\vskip 10.76385pt\cr 0.0pt{\hfil$\displaystyle\qquad K=c_{2}^{2}\sum_{u=0}^{\infty}\tilde{M}^{u/\epsilon}MP_{+}\bigl{\{}\tilde{M}^{u/\epsilon}M\bigr{\}}^{T},\qquad P_{+}=\pmatrix{1&0\cr 0&0\cr},\hfill\hbox to0.0pt{\hss{\tenrm(4.14)}}\crcr}}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\tilde{M}K\tilde{M}^{T}=K-c_{2}^{2}MP_{+}M^{T},\hfill\hbox to0.0pt{\hss{\tenrm(4.15)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad K=\pmatrix{1&0\cr 0&\hat{\Delta}^{-1}}\hfill\hbox to0.0pt{\hss{\tenrm(4.16)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\hat{F}_{r}(\hat{q}_{0},\ldots,\hat{q}_{r})=\Delta\hat{q}_{r}+\hat{F}^{\prime}_{r}(\hat{q}_{0},\ldots,\hat{q}_{r-1}),\hfill\hbox to0.0pt{\hss{\tenrm(4.17)}}\crcr}}$$

$$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr 0.0pt{\hfil$\displaystyle\qquad\pmatrix{\hat{p}_{r}(t)\cr\hat{q}_{r}(t)}=\tilde{M}^{t/\epsilon}\pmatrix{\hat{p}_{r}(0)\cr\hat{q}_{r}(0)}+\sum_{u=\epsilon}^{t}\tilde{M}^{(t-u)/\epsilon}\hat{V}_{r}(\hat{p}_{0}(u),\ldots,\hat{q}_{r-1}(u))\hfill\hbox to0.0pt{\hss{\tenrm(4.18)}}\crcr}}$$

Full text

CERN-TH-2017-057

SMD-based numerical stochastic perturbation theory

Mattia Dalla Bridaa and Martin Lüscherb,c

aDipartimento di Fisica, Università di Milano-Bicocca and

INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

bCERN, Theoretical Physics Department, 1211 Geneva 23, Switzerland

cAlbert Einstein Center for Fundamental Physics

Institute for Theoretical Physics, Sidlerstrasse 5, 3012 Bern, Switzerland

Abstract

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ödinger 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.

1. Introduction

Numerical stochastic perturbation theory (NSPT) [?–?] 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 [?] 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 [?] has recently been briefly looked at in ref. [?] 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.

2. 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},\ldots,q_{n})$ of real coordinates and an action $S(q)$ is considered in this and the following two sections.

2.1 Preliminaries

The action $S(q)$ is assumed to be differentiable and to have an expansion in powers of a coupling $g$ of the form

[TABLE]

where $S_{r}(q)$ is a polynomial in $q$ of degree $d_{r}\geq 2$. Moreover, it is taken for granted that the leading-order term

[TABLE]

is a strictly positive quadratic form in the coordinates.

The observables ${\cal O}(q)$ of interest are assumed to be similarly expandable and their expectation values

[TABLE]

then have a well-defined perturbation expansion with coefficients given by Feynman diagrams as usual.

2.2 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},\ldots,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

[TABLE]

where the momentum $\upsilon$ is randomly chosen from a Gaussian distribution with mean zero and unit variance. The coefficients

[TABLE]

depend on the simulation time step $\epsilon>0$ and a parameter $\gamma>0$ that controls the rotation angle.

In the second step, the molecular-dynamics equations

[TABLE]

are integrated from the current simulation time $t$ to $t+\epsilon$ using a reversible symplectic integration scheme (see subsect. 2.3). The algorithm (momentum rotation followed by the molecular-dynamics evolution) simulates the canonical distribution

[TABLE]

provided the integration errors are negligible or an acceptance-rejection step is included in the update cycle which corrects for these [?]. 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.

2.3 Integration schemes

The molecular-dynamics equations may be integrated by applying a sequence of the elementary steps

[TABLE]

to the current momenta and coordinates, with step sizes $h$ proportional to $\epsilon$. A well-known example is the leapfrog integrator $I_{p,\epsilon/2}I_{q,\epsilon}I_{p,\epsilon/2}$, and several highly efficient schemes are described in ref. [?].

Integrators $I_{\epsilon}$ of this kind are symplectic and they can be (and are here) required to be reversible, i.e. to be such that

[TABLE]

where $P$ stands for the momentum reflection $p\to-p$.

3. Stochastic perturbation theory

Stochastic perturbation theory [?,?] 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.

3.1 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 acceptance-rejection 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

[TABLE]

where the leading-order histories $\hat{p}_{0}(t),\hat{q}_{0}(t)$ coincide with the ones generated by the algorithm in the free theory with action (2.2).

In terms of the coefficients $\hat{p}_{r},\hat{q}_{r}$, the momentum rotation (2.4) becomes

[TABLE]

and the molecular-dynamics integration steps (2.8),(2.9) assume the form

[TABLE]

The forces $\hat{F}_{r}$ in eq. (3.3) are given by

[TABLE]

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 histories $\hat{p}_{r}(t),\hat{q}_{r}(t)$.

3.2 Perturbation expansion of observables

Similarly to the gradient of the action, any observable

[TABLE]

may be expanded in powers of the coupling. The coefficients $k_{r}({\cal O})$ in the perturbation expansion

[TABLE]

of its expectation value (2.3) then coincide with the averages of $\hat{{\cal O}}_{r}(\hat{q}_{0}(t),\ldots,\hat{q}_{r}(t))$ over the simulation time $t$ up to statistical (and integration) errors.

4. 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 $\epsilon$ 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 $\Delta$ in the leading-order part (2.2) of the action.

4.1 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+\epsilon$ amounts to a linear transformation

[TABLE]

of the current momenta and coordinates. The $2n\times 2n$ matrix $M$ in this equation has a block structure,

[TABLE]

with $n\times n$ blocks that are polynomials in $\epsilon$ and the matrix $\Delta$ with some numerical coefficients. In particular, they are commuting real symmetric matrices. In the case of the leapfrog integrator,

[TABLE]

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

[TABLE]

It follows from these relations that the Hamilton function

[TABLE]

is exactly conserved by the integrator. In the following, $\epsilon$ and $\Delta$ are assumed to be such that $M_{qp}$ is non-singular and $\hat{\Delta}$ positive definite. Both conditions are met in the case of the leapfrog integrator if $\epsilon^{2}\left\|\Delta\right\|<4$, where $\left\|\Delta\right\|$ is the largest eigenvalue of $\Delta$. The probability distribution

[TABLE]

is then well-behaved and preserved by the SMD algorithm, since it is preserved both by the momentum rotation and the molecular-dynamics evolution.

4.2 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

[TABLE]

where $\upsilon(u)$ is the momentum chosen randomly in the momentum rotation (2.4) at simulation time $u=0,\epsilon,2\epsilon,\ldots$ and

[TABLE]

The long-time behaviour of the momenta and coordinates thus depends on the properties of the matrix $\tilde{M}$.

The blocks $M_{pp},M_{pq},M_{qp}$ and $\hat{\Delta}$ are commuting real symmetric matrices. Since

[TABLE]

and since $\hat{\Delta}$ is assumed to be positive definite, the eigenvalues of $M_{pp}$ have magnitude strictly less than $1$. The eigenvalues of $\tilde{M}$ are then

[TABLE]

where $\mu$ runs through the eigenvalues of $M_{pp}$. In particular, $|\lambda_{\pm}|<1$ and $\tilde{M}$ 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

[TABLE]

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

[TABLE]

i.e. an inhomogeneous linear equation that has a unique solution, since $\tilde{M}$ is a contraction matrix. A few lines of algebra then show that

[TABLE]

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.

4.3 Convergence beyond the leading order

The assumed structure of the action $S(q)$ implies that the force in the molecular-dynamics integration step (3.3) is of the form

[TABLE]

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

[TABLE]

for all $r\geq 1$. The “vertices” $\hat{V}_{r}(\hat{p}_{0},\hat{q}_{0};\ldots;\hat{p}_{r-1},\hat{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 that $\tilde{M}$ 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 [math] 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 large-time 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.

4.4 Summary

The discussion in this section shows that the SMD algorithm converges to all orders of perturbation theory if the matrix $\Delta$ is strictly positive and if $\epsilon^{2}\left\|\Delta\right\|<\kappa$, where $\kappa$ depends on the molecular-dynamics integrator. In the case of the leapfrog, the 2nd-order OMF and the 4th-order OMF integrators, $\kappa$ is equal to $4$, $6.51$ and $9.87$, respectively (cf. appendix A).

5. 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 ${\tenrm 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.

5.1 Lattice fields

The lattice theory is set up on a $T\times L^{3}$ lattice with periodic boundary conditions in the space directions and Schrödinger functional (SF) [?,?] or open-SF [?] boundary conditions in the time direction. In both cases the link variables $U(x,\mu)$ satisfy

[TABLE]

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. [?]. In particular, all dimensionful quantities are expressed in units of the lattice spacing.

The momenta $\pi(x,\mu)$ of the link variables $U(x,\mu)$ 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,

[TABLE]

includes a conventional weight factor

[TABLE]

which will reappear below in various expressions.

Infinitesimal gauge transformations are fields $\omega(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

[TABLE]

where $\hat{\omega}(x)$ vanishes at $x_{0}=0,T$ and is otherwise unconstrained. A possible choice of scalar product is then

[TABLE]

while in the case of open-SF boundary conditions

[TABLE]

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.

5.2 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

[TABLE]

In particular, the random field in the momentum rotation

[TABLE]

must be distributed with probability density proportional to ${\tenrm e}^{-{1\over 2}(\upsilon,\upsilon)}$.

If the weight factor $w_{x,\mu}$ 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

[TABLE]

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

[TABLE]

Since it will be omitted in perturbation theory, the acceptance-rejection step needed to correct for the integration errors is not described here.

5.3 Perturbation expansion

In perturbation theory, the link variables are represented by a gauge potential $A_{\mu}(x)$ through

[TABLE]

The gauge potential takes values in ${su}(N)$ and satisfies the same boundary conditions as the momentum field

[TABLE]

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 $\hat{\pi}_{0}$ in place of $\pi$) and the integration steps (5.10),(5.11) become

[TABLE]

where $\Delta$ is the symmetric linear operator in the leading-order expression

[TABLE]

for the gauge action.

5.4 Damping of the gauge modes

The space ${\cal H}_{1}$ of gauge potentials may be split into the subspace ${\cal H}_{1}^{L}$ of gauge modes and its orthogonal complement ${\cal H}_{1}^{T}$. 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

[TABLE]

which maps any field $\nu$ in the space ${\cal H}_{0}$ of infinitesimal gauge transformations to a field $d\nu\in{\cal H}_{1}^{L}$. The adjoint operator $d^{\ast}$ goes in the opposite direction and is defined by the requirement that

[TABLE]

for all $\nu\in{\cal H}_{0}$ and $A\in{\cal H}_{1}$. In particular, the subspace ${\cal H}_{1}^{T}$ coincides with the space of gauge potentials $A$ satisfying $d^{\ast}\!A=0$ ($d^{\ast}$ is given explicitly in appendix B).

Since the gauge modes are annihilated by $\Delta$, the ${\cal H}_{1}^{L}$-component of the leading-order gauge field $\hat{U}_{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

[TABLE]

to the fields after each SMD update cycle, where $\omega\in{\cal H}_{0}$ is a new field that is evolved together with the other fields. There are two update steps for this field, the first,

[TABLE]

being applied together with the momentum rotation and the second,

[TABLE]

at the end of the molecular-dynamics evolution. The parameter $\lambda_{0}>0$ controls the feedback from the current gauge field to the gauge-damping field $\omega(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††$\dagger$ In the continuous-time limit $\epsilon\to 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)..

Like the other fields, the gauge-damping field is expanded in a series

[TABLE]

in perturbation theory. At leading order, the new update steps are then

[TABLE]

and the higher-order rules have the familiar hierarchical structure.

5.5 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 matrix $\tilde{M}$ describing the evolution of the fields $(\hat{\pi}_{0},\hat{U}_{0},\hat{\omega}_{0})$ at leading order is a contraction matrix.

Since the gauge modes are zero modes of $\Delta$, the subspace of field vectors

[TABLE]

is invariant under the action of $\tilde{M}$. When restricted to its orthogonal complement (which is invariant too), the matrix describes the evolution of the ${\cal H}_{1}^{T}$-components of $\hat{\pi}_{0},\hat{U}_{0}$. The chosen boundary conditions (SF or open-SF) imply the strict positivity of $\Delta$ in this subspace and $\tilde{M}$ is therefore a contraction matrix there if

[TABLE]

where $\kappa$ 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 of $\tilde{M}$ amounts to applying the linear transformation

[TABLE]

to the fields. The associated eigenvalues are equal to $c_{1}$ or to

[TABLE]

where $\mu^{2}$ is any eigenvalue of $d^{\ast}\kern-1.0ptd$. This operator has no zero modes and some simple estimates then show that $\tilde{M}$ is a contraction matrix in the subspace of gauge modes, provided the bound

[TABLE]

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,

[TABLE]

where $k$ is equal to $1,5/3$ and $3.648$ in the case of the Wilson [?], tree-level Symanzik improved [?,?] and Iwasaki [?] 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 $\epsilon$.

5.6 Inclusion of the quark fields

As in the case of the HMC algorithm [?], 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 $\phi$, of the form

[TABLE]

with $R(U)$ being some gauge-covariant linear operator (see ref. [?], 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††$\dagger$ There is no reason other than simplicity to set the parameter $\gamma$ that determines the rotation angle to the same value for all fields.

[TABLE]

at the beginning of the cycle, where $\eta$ 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 pseudo-fermion 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 leading-order 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.

6. 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 [?,?]. 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 [?] (see ref. [?] for a review).

In the following, the perturbation expansion of the gradient-flow coupling is determined in the ${\tenrm 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 $\overline{\tenrm MS\kern-0.5pt}\kern 0.5pt$ scheme of dimensional regularization is known in infinite volume since a while [?], but a huge effort plus the best currently available techniques were required to be able to extend this calculation to the next order [?]. In finite volume with SF boundary conditions, these techniques do not apply and a similar analytical computation may be practically infeasible at present.

6.1 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 [?]

[TABLE]

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 $\bar{g}^{2}$ coincides with $g_{0}^{2}$ 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\times L$.

On the lattice, the gradient flow is implemented as in ref. [?]. 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^{\tenrm ct}$, for example, denotes the clover expression for the time-like part of the action density).

6.2 Expansion in powers of $\alpha_{s}$

The gradient flow coupling has an expansion

[TABLE]

in powers of the running coupling $\alpha_{s}$ in the $\overline{\tenrm MS\kern-0.5pt}\kern 0.5pt$ scheme of dimensional regularization at momentum scale $q$, with coefficients $k_{1},k_{2},\ldots$ depending on $q$ and $L$. If $q$ is scaled with $L$ like

[TABLE]

the coefficients are of order $1$ and independent of $L$ in the continuum theory [?].

In the following $k_{1}$ and $k_{2}$ are computed by combining the expansion

[TABLE]

obtained in NSPT with the relation

[TABLE]

between $\alpha_{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 [?]. 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.

6.3 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 [?]††$\dagger$ In ref. [?] the improvement coefficient $c_{\eightrm t}$ is denoted by $c_{\eightrm G}^{\prime}$.

[TABLE]

In the case of the Wilson action,

[TABLE]

as was shown long ago [?–?]. The counterterms lead to further terms in the forces that drive the molecular-dynamics evolution, but do not require a modification of the general framework described in sect. 5. Alternatively, the effects of the counterterms can be included in the calculations by treating $c_{\tenrm 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,\mu)$ at flow time $t$ is represented by a gauge potential $B_{\mu}(t,x)$ according to

[TABLE]

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.

Table 1. Simulation parameters

[TABLE]

6.4 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 $\lambda_{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 $\Delta t_{\tenrm ms}/\epsilon$ update cycles using a 3rd-order Runge–Kutta integrator for the gradient flow [?], with step size varying from $0.002$ at small flow times to $0.1$ at large times. With this scheme, the gradient-flow integration errors are guaranteed to be completely negligible with respect to the statistical errors. The number $N_{\tenrm 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\in\{0.2,0.3,0.4\}$ and all choices w,c,ws,$\ldots$ of the action density (cf. subsect. 6.1).

7. 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 $\epsilon$ and the SMD parameter $\gamma$. 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. [?,?] for related complementary studies of the $\phi^{4}$ theory).

7.1 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 $\hat{{\cal O}}_{r}$ of an observable ${\cal O}$ in general does not coincide with the order-$2r$ coefficient of another observable. The time average of $\hat{{\cal O}}_{r}^{2}$ and thus the variance of $\hat{{\cal O}}_{r}$ are then not necessarily determined by the theory alone. For illustration, the dependence on $\gamma$ of the variances of the coefficients in

[TABLE]

(where $x_{0}=L/2$ and $\sqrt{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 $\gamma$, which is therefore entirely an effect of the change of algorithm.

Fig. 1. Variances of $\hat{E}^{\ninerm c}_{k}$ at fixed $L=16$, $c_{\ninerm t}=1$, $c=0.3$ and $\epsilon=0.238$ versus $\gamma$. The dotted lines connecting the data points are drawn to guide the eye. Beyond $\gamma=3$ the variances are practically constant.

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

[TABLE]

of the integrated autocorrelation times and variances of the order-$r$ coefficients $\hat{{\cal O}}_{r}$ of the observables ${\cal 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 $\hat{E}_{k}^{\tenrm c}$, for example, grow monotonically with $\gamma$ (see table 2). At the chosen point in parameter space, the associated products (7.2) are then minimized at values of $\gamma$ 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 $\gamma$, but a range of values, where all coefficients of interest are obtained reasonably efficiently.

Table 2. Autocorrelation times of $\hat{E}_{k}^{\ninerm c}$ and associated products (7.2)$\dagger$

[TABLE]

7.2 Critical slowing down

The behaviour of the autocorrelation times and variances near the continuum limit $L\to\infty$ 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}$ [?,?], while their variances grow at most logarithmically [?].

At large $\gamma$, the variant of NSPT studied here is expected to behave similarly, since the SMD algorithm then effectively integrates the Langevin equation. Choosing $\gamma$ to 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 $\gamma$ is scaled like $1/L$ [?]. 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 $\hat{E}_{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,\ldots,40$ in table 1) were chosen so that subsequent measurements are practically uncorrelated. At fixed $\gamma=3$ the required increase with $L$ of the measurement time separation $\Delta t_{\tenrm ms}$ then turned out to be quite moderate. Moreover, from $L=24$ to $L=40$, the variances of $\hat{E}_{k}$ grow only slowly (at $c=0.3$ and for $k=0$, $1$ and $2$ by about $2$, $19$ and $30$ percent).

7.3 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 $\hat{E}_{0}$ calculated in the runs listed in table 1 in fact all agree with the known analytic formula [?,?] within a relative statistical precision of about $2\times 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 $\hat{H}_{0},\ldots,\hat{H}_{4}$ of the Hamilton function $H$ are accurately conserved, with deficits decreasing like $\epsilon^{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\geq 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.

Fig. 2. Dependence of $k^{\ninerm c}_{1}$ and $k^{\ninerm c}_{2}$ at $L=24$, $c_{\ninerm t}=1$ and $c=0.3$ on the simulation step size $\epsilon$ (note the log scale on the abscissa). In the range of data shown, $\epsilon$ increases from $0.1$ to $0.238$. Each data point is based on approximately $10^{4}$ statistically independent measurements of the coefficients $\hat{E}^{\ninerm c}_{0}$, $\hat{E}^{\ninerm c}_{1}$ and $\hat{E}^{\ninerm c}_{2}$.

7.4 Extrapolation to the continuum limit

With ${\tenrm O}(a)$-improvement in place, the $L$-dependence of the one-loop coefficient $k_{1}^{\tenrm w}$ is asymptotically given by [?,?]

[TABLE]

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\geq 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_{1}^{\tenrm c},k_{1}^{\tenrm ws},\ldots$ look much the same.

Fig. 3. Extrapolation of the simulation results for $k_{1}^{\ninerm w}$ at $c=0.3$ to the continuum limit (open point). The full line is obtained by fitting the data in the range $10\leq L\leq 40$ with the asymptotic expression (7.3). A linear fit in the range $16\leq L\leq 40$ (dashed grey line) lies practically on top of the full line.

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 standard 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.

Fig. 4. Extrapolation of the simulation results for $k_{2}^{\ninerm w}$ at $c=0.3$ to the continuum limit (open point). The full line is obtained by fitting the data in the range $10\leq L\leq 40$ according to eq. (7.4), with the two lowest modes of the design matrix projected away. In the range $16\leq L\leq 40$ the data can also be fitted by a straight line (dashed grey line).

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

[TABLE]

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.

7.5 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.

Table 3. Values of $k_{1}$ and $k_{2}$ in the continuum limit

[TABLE]

Infinite-volume limit. The gradient-flow coupling in infinite volume runs with the flow time $t$ and may be expanded in powers of $\alpha_{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)$ [?] and $k_{2}=-0.9822(1)$ [?]. 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},\ldots,k^{\tenrm t}_{2}$ at $c=0.2$ are actually already quite close to the infinite-volume values quoted above.

Three-loop $\beta$-function. The $L$-dependence of the gradient-flow coupling $\alpha=\bar{g}^{2}/4\pi$ is governed by the renormalization group equation

[TABLE]

where $\beta_{0}=11(2\pi)^{-1}$ and $\beta_{1}=51(2\pi)^{-2}$ are the universal one- and two-loop coefficients of the Callan–Symanzik $\beta$-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 $\overline{\tenrm MS\kern-0.5pt}\kern 0.5pt$ 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.

8. 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 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.

Table 4. Ratio $\rho_{2}=\beta_{2}/\beta_{0}$ of coefficients of the $\beta$-function

[TABLE]

The SMD algorithm provides a technically attractive basis for NSPT. Compared to the traditional version of NSPT [?–?], 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 molecular-dynamics equations allow the integration errors to be easily kept under control. For the reasons mentioned above, some tuning of the friction parameter $\gamma$ 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 [?]. 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 $\alpha_{s}^{l+1}$ in the perturbation series (6.2). An extension of the computations to three-loop 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 $\alpha_{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.

Appendix A. OMF molecular-dynamics integrators

Among the popular integrators used in lattice QCD, there are two schemes proposed by Omelyan, Mryglod and Folk (OMF) [?], one of second order in the integration step size and the other of fourth order. Here their efficiency is studied by comparing the distribution (4.8) actually simulated at leading order of stochastic perturbation theory with the exact distribution.

A.1 2nd-order OMF integrator

The sequence of update steps (2.8),(2.9) is

[TABLE]

in this case. Using an algebraic manipulation program, the expressions

[TABLE]

are then easily obtained. The convergence of the SMD algorithm is thus guaranteed if $\epsilon^{2}\left\|\Delta\right\|<6.51$. With respect to the leapfrog integrator, and if the step size $\epsilon$ is adjusted so that the number of force evaluations per unit time is the same, this scheme achieves a reduction of the relative deviation of $\hat{\Delta}$ from $\Delta$ by a factor $25$.

A.2 4th-order OMF integrator

This integrator is of the form (A.1) too, but with $11$ steps instead of $5$. The associated step sizes in units of $\epsilon$ are

[TABLE]

and the matrices $M_{pq}$ and $M_{qp}$ are polynomials in $\epsilon^{2}\Delta$ of degree $5$ and $4$ times $\epsilon\Delta$ and $\epsilon$, respectively. In this case, $M_{qp}$ is non-singular and

[TABLE]

is positive if $\epsilon^{2}\left\|\Delta\right\|<9.87$.

As is evident from eqs. (A.6),(A.7), this integrator is impressively accurate. If $\left\|\Delta\right\|=16$ and $\epsilon=0.2$, for example, the relative deviation of $\hat{\Delta}$ from $\Delta$ is at most $1.6\times 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 $\epsilon=0.08$).

Appendix B. Explicit form of $\tenmib d^{\ast}$

The operator $d^{\ast}$ depends on the chosen boundary conditions and the scalar products in the field spaces ${\cal H}_{0}$ and ${\cal H}_{1}$ (cf. subsect. 5.1). For SF boundary conditions

[TABLE]

where

[TABLE]

while in the case of open-SF boundary conditions,

[TABLE]

where now

[TABLE]

Appendix C. Gauge-damped stochastic equations

In the continuous-time limit $\epsilon\to 0$, the gauge-damped SMD algorithm described in sect. 5 integrates the stochastic equations

[TABLE]

where $\pi_{t},U_{t},\omega_{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

[TABLE]

and

[TABLE]

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 $\pi_{t},U_{t}$ and $\eta_{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 [?]. The gauge damping can also easily be shown to coincide with the standard one in the Langevin limit $\gamma\to\infty$ [?].

Appendix D. Tables of simulation results

Table 5. Values of the expansion coefficient $k_{1}$ at $c=0.2$

[TABLE]

Table 6. Values of the expansion coefficient $k_{1}$ at $c=0.3$

[TABLE]

Table 7. Values of the expansion coefficient $k_{1}$ at $c=0.4$

[TABLE]

Table 8. Values of the expansion coefficient $k_{2}$ at $c=0.2$

[TABLE]

Table 9. Values of the expansion coefficient $k_{2}$ at $c=0.3$

[TABLE]

Table 10. Values of the expansion coefficient $k_{2}$ at $c=0.4$

[TABLE]

References

[1] F. Di Renzo, G. Marchesini, P. Marenzoni, E. Onofri, Lattice perturbation theory on the computer, Nucl. Phys. Proc. Suppl. 34, 795 (1994).

[2] F. Di Renzo, E. Onofri, G. Marchesini, P. Marenzoni, Four loop result in SU(3) lattice gauge theory by a stochastic method: Lattice correction to the condensate, Nucl. Phys. B426, 675 (1994).

[3] F. Di Renzo, L. Scorzato, Numerical stochastic perturbation theory for full QCD, JHEP 0410 (2004) 073.

[4] S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Hybrid Monte Carlo, Phys. Lett. B195 (1987) 216.

[5] A. M. Horowitz, Stochastic quantization in phase space, Phys. Lett. 156B (1985) 89; The second order Langevin equation and numerical simulations, Nucl. Phys. B280 [FS18] (1987) 510; A generalized guided Monte Carlo algorithm, Phys. Lett. B268 (1991) 247.

[6] M. Dalla Brida, M. Garofalo, A. D. Kennedy, Numerical stochastic perturbation theory and gradient flow in $\phi^{4}$ theory, PoS LATTICE2015 (2016) 309.

[7] K. Jansen, C. Liu, Kramers equation algorithm for simulations of QCD with two flavors of Wilson fermions and gauge group SU(2), Nucl. Phys. B453 (1995) 375 [Erratum: ibid. B459 (1996) 437].

[8] I. P. Omelyan, I. M. Mryglod, R. Folk, Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations, Comp. Phys. Commun. 151 (2003) 272.

[9] G. Parisi, Y.-S. Wu, Perturbation theory without gauge fixing, Sci. Sin. 24 (1981) 483.

[10] P. H. Damgaard, H. Hüffel, Stochastic quantization, Phys. Rept. 152 (1987) 227.

[11] M. Lüscher, R. Narayanan, P. Weisz, U. Wolff, The Schrödinger functional --- a renormalizable probe for non-Abelian gauge theories, Nucl. Phys. B384 (1992) 168.

[12] S. Sint, On the Schrödinger functional in QCD, Nucl. Phys. B421 (1994) 135.

[13] M. Lüscher, Step scaling and the Yang--Mills gradient flow, JHEP 1406 (2014) 105.

[14] K. G. Wilson, Confinement of quarks, Phys. Rev. D10 (1974) 2445.

[15] P. Weisz, Continuum limit improved lattice action for pure Yang--Mills theory (I), Nucl. Phys. B212 (1983) 1.

[16] M. Lüscher, P. Weisz, On-shell improved lattice gauge theories, Commun. Math. Phys. 97 (1985) 59 [Erratum: ibid. 98 (1985) 433].

[17] Y. Iwasaki, Renormalization group analysis of lattice theories and improved lattice action. II -- Four-dimensional non-Abelian SU(N) gauge model, preprint UTHEP-118 (1983) [arXiv:1111.7054v1].

[18] M. Lüscher, Computational strategies in lattice QCD, in: Modern perspectives in lattice QCD, eds. L. Lellouch et al. (Oxford University Press, New York, 2011).

[19] P. Fritzsch, A. Ramos, The gradient flow coupling in the Schrödinger functional, JHEP 1310 (2013) 008.

[20] M. Dalla Brida et al. (ALPHA Collab.), Slow running of the gradient flow coupling from 200 MeV to 4 GeV in $N_{f}=3$ QCD, Phys. Rev. D95 (2017) 014507.

[21] I. Campos et al., Non-perturbative running of quark masses in three-flavour QCD, PoS LATTICE2016 (2016) 201.

[22] M. Lüscher, P. Weisz, U. Wolff, A numerical method to compute the running coupling in asymptotically free theories, Nucl. Phys. B359 (1991) 221.

[23] M. Lüscher, Advanced lattice QCD, in: Probing the Standard Model of Particle Interactions (Les Houches 1997), eds. R. Gupta et al. (Elsevier, Amsterdam, 1999).

[24] M. Lüscher, Properties and uses of the Wilson flow in lattice QCD, JHEP 1008 (2010) 071 [Erratum: ibid. 1403 (2014) 092].

[25] R. V. Harlander, T. Neumann, The perturbative QCD gradient flow to three loops, JHEP 1606 (2016) 161.

[26] M. Lüscher, P. Weisz, Perturbative analysis of the gradient flow in non-Abelian gauge theories, JHEP 1102 (2011) 051.

[27] M. Lüscher, P. Weisz, Two-loop relation between the bare lattice coupling and the $\overline{\tenrm MS\kern-0.5pt}\kern 0.5pt$ coupling in pure ${\tenrm SU}(N)$ gauge theories, Phys. Lett. B349 (1995) 165.

[28] M. Lüscher, R. Sommer, P. Weisz, U. Wolff, A precise determination of the running coupling in the SU(3) Yang--Mills theory, Nucl. Phys. B413 (1994) 481.

[29] A. Bode, U. Wolff, P. Weisz, Two-loop computation of the Schrödinger functional in pure SU(3) lattice gauge theory, Nucl.Phys. B540 (1999) 491.

[30] A. Bode, P. Weisz, U. Wolff, Two-loop computation of the Schrödinger functional in lattice QCD, Nucl. Phys. B576 (2000) 517 [Errata: ibid. B600 (2001) 453, ibid. B608 (2001) 481].

[31] M. Dalla Brida, M. Garofalo, A. D. Kennedy, An investigation of new methods for numerical stochastic perturbation theory in $\varphi^{4}$ theory, arXiv:1703.04406 [hep-lat]

[32] J. Zinn–Justin, Renormalization and stochastic quantization, Nucl. Phys. B275 [FS17] (1986) 135.

[33] J. Zinn–Justin, D. Zwanziger, Ward identities for the stochastic quantization of gauge fields, Nucl. Phys. B295 [FS21] (1988) 297.

[34] M. Lüscher, Statistical errors in stochastic perturbation theory, notes (2015), http://luscher.web.cern.ch/luscher/notes/enspt.pdf.

[35] K. Symanzik, Cutoff dependence in lattice $\phi^{4}_{4}$ theory, in: Recent developments in gauge theories (Cargèse 1979), eds. G. ’t Hooft et al. (Plenum, New York, 1980).

[36] K. Symanzik, Some topics in quantum field theory, in: Mathematical problems in theoretical physics, eds. R. Schrader et al., Lecture Notes in Physics, Vol. 153 (Springer, New York, 1982).

[37] D. Zwanziger, Covariant quantization of gauge fields without Gribov ambiguity, Nucl. Phys. B192 (1981) 259.