# Perturbative expansions from Monte Carlo simulations at weak coupling: Wilson loops and the static-quark self-energy

###### Abstract

Perturbative coefficients for Wilson loops and the static-quark self-energy are extracted from Monte Carlo simulations at weak coupling. The lattice volumes and couplings are chosen to ensure that the lattice momenta are all perturbative. Twisted boundary conditions are used to eliminate the effects of lattice zero modes and to suppress nonperturbative finite-volume effects due to phases. Simulations of the Wilson gluon action are done with both periodic and twisted boundary conditions, and over a wide range of lattice volumes (from to ) and couplings (from to ). A high precision comparison is made between the simulation data and results from finite-volume lattice perturbation theory. The Monte Carlo results are shown to be in excellent agreement with perturbation theory through second order. New results for third-order coefficients for a number of Wilson loops and the static-quark self-energy are reported.

## I Introduction

Simulations using highly-improved lattice actions have become commonplace in recent years. Effective use of these requires perturbative matching calculations for masses, coupling constants and currents, among other quantities. Higher-order perturbative calculations for these actions are laborious but they are essential in order to obtain precision results for most observables.

An alternative to doing calculations in analytical perturbation theory is to directly measure short-distance quantities in Monte Carlo simulations at weak coupling, as proposed in Refs. [1, 2]. One exploits the fact that the lattice theory on a finite volume enters a perturbative phase at weak coupling. In effect the couplings and volumes in the simulations are chosen to ensure that the lattice momenta are all perturbative (up to possible zero modes). In this way one can in principle extract perturbative expansions for many quantities, by fitting Monte Carlo data for appropriate correlation functions to a power series in the coupling.

This approach has been shown to reproduce analytical results for the first-order mass renormalization for Wilson fermions, and the first-order additive self-energy for NRQCD fermions [1, 2]. In addition, preliminary estimates of some third-order Wilson loop coefficients were made in Ref. [3]. An extension of this technique to background field calculations was considered in Ref. [4]. Preliminary work on perturbative simulations of quark actions has also been done [5, 6].

The Monte Carlo method as implemented here requires as input from conventional lattice perturbation theory only the expansion of the average plaquette, and of the static potential (or some other quantity that defines a physical coupling constant), to the desired order, along with an estimate of the scale relevant to the quantity of interest. These inputs are necessary because we use a renormalization scheme defined through the perturbative expansion of the static-quark potential [7, 8]. The renormalized coupling in this scheme can be extracted from measured values of the average plaquette, given its perturbative expansion.

The method proceeds as follows. Simulations are done at a several different values of ( for SU() gauge theory). At each we use the measured value of the plaquette to solve for the value of the renormalized coupling , at the scale that is optimal for the plaquette [7]. We then run the couplings at each to the scale appropriate to the quantity of interest, whose expectation values are then fit to a truncated series in . The fit yields numerical values for the perturbative coefficients. To assess the effects of the truncation of the perturbation series at finite order in , the fits are done including many higher-order terms, beyond the order of interest, but where the fits incorporate constraints on the coefficients [9], which are required to lie in a range of values that is consistent with a well-behaved perturbative expansion. One can also improve the quality of the results by using lower-order coefficients from conventional perturbation theory, if available, in order to further constrain the fits to the Monte Carlo data, thereby obtaining more accurate values for previously unknown higher-order terms.

In order to ensure that the lattice momenta are all sufficiently perturbative, simulations are done on lattices for which

(1) |

where is the QCD scale parameter, is the lattice spacing, and is the physical length of the lattice. In order to minimize perturbative finite-volume errors we must also ensure that the spacing between lattice momenta is small compared to the characteristic momentum scale associated with the quantity of interest

(2) |

In practical simulations the lattices are such that Eq. (1) is extremely well satisfied: ranges from to in the analyses in this paper. Moreover the quantities studied here all have characteristic scales near the ultraviolet cutoff, hence is independent of , and Eq. (2) can be easily satisfied. In theories with additional scales, such as quark masses , the characteristic scale is proportional to , where is some dimensionless function, hence is independent of provided that is adjusted to keep fixed. Perturbative finite-volume effects can also be analyzed in detail by running simulations at several different volumes.

In simulations with periodic boundary conditions (PBC) however there are lattice zero modes which will violate Eq. (1). One must also ensure that nonperturbative finite-volume effects arising from phases are sufficiently suppressed. One way to achieve this, using a local updating algorithm, is to work on lattices with sufficiently large volumes, where the simulation is started with perturbative initial values for the links (e.g., by setting all links to the identity). A more powerful approach, analyzed in detail here, is to adopt boundary conditions that eliminate zero modes and suppress transitions between different phases.

It is also desirable to have a more general means for dealing with potential infrared problems. The effects of lattice zero modes on most observables is not known and, while these could in some cases be accounted for by a numerical fit to the data, it is preferable to eliminate these states from the outset. More generally zero modes can significantly alter the expected perturbative form of correlation functions. For example, the presence of zero modes would pose a significant problem for extracting quark masses, as one could not assume the existence of an isolated pole in the perturbative quark propagator. Twisted boundary conditions (TBC) [10, 11, 12] can be used to eliminate zero modes, and are easily incorporated into simulations using existing code for a given action.

In this work we present a comprehensive study of this Monte Carlo method for extracting perturbative quantities [13]. We do simulations of SU(3) gauge theory using the Wilson gluon action. The evolution of the coupling is known for this theory through three-loop order which, in principle, allows one to use the method to determine perturbative expansions through fourth order. Simulations are done for both PBC and TBC. We show that using TBC is an effective means of suppressing nonperturbative finite-volume effects due to phases, as well as eliminating the effects of lattice zero modes. We also make an extensive analysis of perturbative finite-volume effects for TBC. We analyze simulation results for a large set of Wilson loops and the static-quark self-energy, over a wide range of couplings and lattice volumes. A high-precision comparison is made between the simulation data and results from finite-volume lattice perturbation theory. The Monte Carlo results are shown to be in excellent agreement with perturbative calculations through second order, which are available for these observables for both periodic [14, 15] and twisted [16] boundary conditions. New results for third-order coefficients for fourteen different Wilson loops and the static-quark self-energy are also reported.

Wilson loops provide a good quantity for a first test of the Monte Carlo method, since small loops are relatively insensitive to nonperturbative phases, and have small finite-volume errors. The perturbative expansion of small Wilson loops is also relevant to determinations of the strong coupling from lattice simulations [8]. The calculation of the static-quark self-energy is considerably more involved, as it is very sensitive to nonperturbative phases and has large perturbative finite-volume corrections. The static-quark self-energy thus represents a good prototype for more realistic calculations of other perturbative quantities, such as quark masses. The self-energy is also useful for determinations of the -quark mass from lattice simulations [17].

The rest of this paper is organized as follows. In Sect. II we use PBC to study Wilson loops at large . Simulations are done on lattices at nine couplings. We evaluate Wilson loops of sizes with . The results are fit to a truncated perturbation series, using the renormalized coupling evaluated at a scale that is appropriate to each Wilson loop [7]. We explicitly subtract the leading effects of lattice zero modes from the Monte Carlo data using an existing analytical calculation [18]. A detailed comparison is made with perturbation theory through second order, and estimates are made of the third-order coefficients.

In Sect. III we review how simulations are done using twisted boundary conditions, and their use in analytical perturbation theory. We demonstrate that these boundary conditions can be used to virtually eliminate nonperturbative effects due to phases even on lattices with very small volumes.

The static-quark self-energy is analyzed in Sect. IV. We extract the self-energy from the gauge-invariant Polyakov line, which describes the propagation of a static quark across the entire time-extent of the lattice. Other extractions of the self-energy have relied upon large Wilson loops [15, 19], and the gauge-fixed quark propagator [5] (the latter methods could prove more efficient in Monte Carlo simulations, as one can limit the propagation time on a lattice of a given size). We make a perturbative analysis of the self-energy on finite lattices, and show that one can sum leading logarithms in the finite-volume correction. We present results of simulations of the self-energy, with runs at nine volumes at each of nine couplings. We compare the simulation results with perturbation theory through second order, over the whole range of lattice sizes studied here, and make an estimate of the third-order self-energy, including an extrapolation to the infinite-volume limit.

During the course of the calculations described in this paper, we also investigated the question of how to design the most efficient calculations through intelligent parameter choices, by using the techniques of constrained curve fitting [9]. In calculations like those in this paper, for example, one can choose the couplings in order to minimize the simulation cost required to achieve a given precision in the perturbative coefficients. Although the Monte Carlo calculations described in this paper were designed before these optimization techniques were worked out, we include a description of them here, in an Appendix, for use in future calculations.

Some conclusions and prospects for future work are briefly discussed in Sect. V.

## Ii Wilson loops

In this section we analyze simulations of Wilson loops at large . The Wilson gauge-field action for SU(3) color is used, where

(3) |

and is the plaquette. Simulations were done on lattices at nine couplings. Details of the simulation parameters are given in Table 1. Periodic boundary conditions were used here, in order to make a direct comparison with the first- and second-order perturbative coefficients calculated on finite lattices by Heller and Karsch [14]. We will also verify that the simulations reproduce the effects of lattice zero modes to leading order.

Measurements | |||||
---|---|---|---|---|---|

0.010 | 47.746 | 2320 | 0.957542(1) | 0.01049 | |

0.015 | 31.831 | 2459 | 0.935857(2) | 0.01614 | |

0.020 | 23.873 | 2112 | 0.913829(3) | 0.02209 | |

0.025 | 19.099 | 1558 | 0.891441(3) | 0.02839 | |

0.030 | 15.915 | 860 | 0.868599(9) | 0.03510 | |

0.035 | 13.642 | 746 | 0.845305(8) | 0.04225 | |

0.040 | 11.937 | 748 | 0.821472(9) | 0.04992 | |

0.045 | 10.610 | 500 | 0.797038(11) | 0.05819 | |

0.050 | 9.459 | 500 | 0.771872(11) | 0.06719 |

We analyze the logarithm of the Wilson loop

(4) |

using a renormalized coupling that is determined from measured values of the plaquette, according to [7, 8]

(5) |

The coupling is defined such that the logarithm of the plaquette has no third- or higher-order terms in its perturbative expansion [8]. Other quantities, of course, do have higher-order terms when expressed as a series in . One can express as a series in the bare lattice coupling , using the third-order expansion of the plaquette given in Ref. [20]

(6) |

The large coefficients in this expansion are an artifact of ; using eliminates large renormalizations of the bare coupling. We also note that the logarithm of the Wilson loop is better behaved perturbatively than the Wilson loop itself, due to the exponentiation of the perturbative perimeter law; this is also why we have divided by the perimeter in defining the perturbative coefficients in Eq. (4).

The perturbation series for each Wilson loop is evaluated using the renormalized coupling at a scale determined according to the procedure of Ref. [7]; the scale corresponds to the typical momentum carried by a gluon in the leading-order diagram for a given quantity. The scales are given in Table 2. The couplings at these scales are evaluated by first measuring the average plaquette in the simulation, and solving Eq. (5) for . We then evolve to the scale appropriate to the quantity under consideration using the universal second-order beta function, plus the third-order term for , with [21]

(7) | |||||

For our quenched simulations , , and , where is the third beta function coefficient in the scheme, and where can be obtained from existing three-loop calculations [20, 22], as described in Ref. [8]. The values of and for our simulations are given in Table 1.

Loop | Loop | ||
---|---|---|---|

3.07 | 2.46 | ||

3.01 | 2.45 | ||

2.96 | 2.38 | ||

2.95 | 2.35 | ||

2.65 | 2.30 | ||

2.56 | 2.27 | ||

2.49 | 2.23 |

One can convert a perturbative expansion in to one in , at the scale appropriate to a particular quantity, through third order, using

(8) | |||||

This connection can be obtained from the third-order expansion of the plaquette in the bare coupling [20], which can be used to solve for in terms of , given its definition in Eq. (5); one can then use a perturbative expansion of the evolution equations [23] to eliminate in favor of . Equation (8) extends the second-order connection between and given in Ref. [7].

The first- and second-order perturbative coefficients for the Wilson loop were computed by Heller and Karsch [14] for an expansion in the bare lattice coupling. We convert this expansion to a series in using Eq. (8). These perturbative calculations were done on finite lattices with periodic boundary conditions, neglecting the contribution of lattice zero modes. Our simulations of the Wilson loops were also done with PBC but do contain the effects of zero modes. In Sect. III we will consider simulations using twisted boundary conditions to eliminate these states. For Wilson loops however we can make use of an analytical calculation of the zero mode piece of the first-order coefficient in Eq. (4), due to Coste et al. [18]

(9) |

where is the lattice volume. We will use this expression to explicitly subtract the leading-order effects of the zero modes from the Monte Carlo data. In the following we will show the Monte Carlo data after first making this zero mode subtraction, unless explicitly noted otherwise; hereafter we will use to denote the (finite-volume) coefficients without zero mode contributions.

We present Monte Carlo results for the Wilson loop in Fig. 1, in terms of the quantity

(10) |

which should exhibit the limit as . We extract estimates of the perturbative coefficients from the Monte Carlo data by fitting the results to the series expansion Eq. (4) where, to begin with, we treat all coefficients as unknown. We will compare the fit values for and with the results from analytical perturbation theory, which provides a stringent test of the Monte Carlo method.

An important aspect of the fitting procedure is how to reliably account for the systematic error arising from the truncation of the fit function at a finite order in . Including too few terms in the expansion in results in a poor fit to the data, while including too many higher-order terms results in very poorly constrained values for the lowest-order coefficients which should, in fact, make the dominant contributions to the data. This situation can be remedied by incorporating constraints on the coefficients, which are required to lie in a range of values that is consistent with our expectation that the perturbative expansion is well behaved. We do this in practice by using conventional least-squares fitting routines, where the is augmented according to:

(11) |

which tends to constrain the fit values for the to the interval . This approach can be motivated by Bayesian statistical analysis [9].

If perturbation theory is reliable we expect the coefficients to be of . We performed least-squares fits to Eq. (4), minimizing with and for the first five orders in the expansion. The dashed line in Fig. 1 shows the results of the fit for the Wilson loop; note that the curvature in the Monte Carlo data shows the sensitivity of the simulations to the third-order term in the perturbative series. The quality of the fits is very good, with -values in excess of 50%.

Loop | |||||
---|---|---|---|---|---|

1.2037(2) | 1.2039 | 0.0(5) | |||

1.2587(2) | 1.2589 | 0.4(5) | |||

1.4337(2) | 1.4338 | 1.1(5) | |||

1.6089(3) | 1.6089 | 2.5(6) | |||

1.7067(4) | 1.7067 | 3.4(6) | |||

1.7693(6) | 1.7690 | 4.3(7) |

The measured values of and are in excellent agreement with perturbation theory, as shown in Table 3, with an accuracy of a few parts in for the first-order coefficients and a few parts in for the second-order coefficients. The third-order coefficient can also be resolved, here with almost no input from analytical perturbation theory. The fit values are very stable to changes in the values of and used in Eq. (11).

Loop | Loop | ||
---|---|---|---|

0.43(9) | 2.52(17) | ||

0.66(11) | 2.53(15) | ||

0.84(12) | 2.98(17) | ||

0.94(14) | 3.26(19) | ||

1.41(11) | 3.40(19) | ||

1.91(13) | 3.71(21) | ||

2.28(15) | 3.91(23) |

Note that if the Monte Carlo data are fit with constrained to its value from perturbation theory, then the errors on are reduced by a factor of about three, with fit values in agreement with perturbation theory within the reduced errors. Similarly, we obtain more accurate results for by fitting the Monte Carlo data with and constrained to their perturbative values. We did fits to Eq. (4) for the next three orders in the perturbative expansion, minimizing using and for . The results for are given in Table 4, where the fit errors are seen to be about 10%.

It is also interesting to verify that the simulations reproduce the leading effects of the zero modes. A convenient way to visualize these effects is to plot the quantity

(12) |

where the first-order coefficient is set to its perturbative value. We plot for the loop in Fig. 2, after subtracting the leading-order zero mode term from the Monte Carlo data. We see that the data reproduce the second-order coefficient from perturbation theory, with as . In Fig. 3 we plot but where the leading-order zero mode is not subtracted from the data. We see evidence of singular behavior in at small coupling, indicating that the first-order term is not completely removed from the Monte Carlo data when the zero mode component is not treated. The dashed line in Fig. 3 shows the results of a fit to Eq. (4), taking account of the leading zero mode contribution, where the term is included in the fit line. This shows explicitly that the Monte Carlo data are sufficiently accurate to reveal the small contribution from the zero modes, at sufficiently small couplings.

We also present results for the residual

(13) |

for the loop in Fig. 4, which is convenient for visualizing the sensitivity of the Monte Carlo data to the third- and fourth-order terms. The statistical errors in the Monte Carlo data are too large to resolve , although the best fits suggest that is of the same order as for all the Wilson loops analyzed here.

A potential complication in our analysis of is that we have only corrected the Wilson loop data for the effects of zero modes to first order. However we expect that the leading contribution from zero modes that remains is of which, given the lattice volume and the range of couplings analyzed here, should only be comparable to terms of ). In fact there is no visible effect of zero modes beyond first order, within statistical errors; this would should up as singular behavior in at small (compare Fig. 4 for with the two plots of in Figs. 2 and 3). We also note that while we have extracted values for on a finite lattice, the volume is large enough that the results should give a good approximation to the coefficients on an infinite lattice (with the corrections expected to be of ).

A determination of higher order terms in the expansion of the and Wilson loops has also been made in Ref. [24], using numerical simulations of the Langevin equations, where a perturbative expansion in the bare lattice coupling is applied to the evolution equations themselves. The results are presented in Ref. [24] as an expansion of the Wilson loops in powers of the bare lattice coupling, . We can convert our third-order result for the expansion of in , to an expansion of in , by using the inverse of Eq. (8). We find for the Wilson loop, in agreement with the result reported in Ref. [24].

We note parenthetically that the expansion of the Wilson loop itself is very poorly convergent, and that the vanishingly small value of for the loop is accidental. One finds very large expansion coefficients for other Wilson loops. These expansions are tamed by taking the logarithm, and expressing the series in . For example, our results give for the loop, with this large value arising almost entirely from the exponentiation of the perturbative perimeter law, and the renormalization of the bare coupling (compare with for the logarithm of the loop). In Sect. IV we show that the third-order expansion of the static-quark self-energy is also very reasonable when expressed in terms of , but is very poorly convergent when expressed in terms of .

## Iii Twisted boundary conditions

### iii.1 Formalism

The analysis of the preceding section shows that simulations at large can be used to make accurate determinations of perturbative quantities at higher orders than have been achieved using conventional perturbation theory. However, as discussed in Sect. I, simulations with periodic boundary conditions (PBC) are subject to the effects of lattice zero modes, and a more convenient method for dealing with potential infrared problems in more general situations is required. Fortunately lattice zero modes can be completely eliminated by using twisted boundary conditions (TBC). We also find that TBC significantly reduce nonperturbative finite-volume effects due to phases.

Twisted boundary conditions [10, 11, 12] for the link fields resemble a gauge transformation on fields which cross through selected lattice boundaries

(14) |

where we take the twist matrices to be constant. A twist must be applied to at least two boundaries and , as the transformation matrix for a single boundary can be completely eliminated by a field redefinition. The requirement that , which can be connected to by crossing the two lattice boundaries in different orders, be single-valued implies that the twist matrices must satisfy the algebra

(15) |

where we consider the general SU() gauge theory in most of this section. A pair of twist matrices generates a multiplication table that forms a discrete subgroup of SU(). In particular [11]

(16) |

where is the unit matrix.

The Wilson gluon action for TBC is written in terms of link variables in the usual way, Eq. (3). The matrices that generate gauge transformations

(17) |

are subject to the same TBC as the links (Eq. (14)). One consequence of this is that the Polyakov line in a twisted direction must have an additional factor of the corresponding twist matrix, if it is to be gauge-invariant:

(18) |

The only zero-action fields are pure-gauge configurations [11], including possible phases. The action also possesses a discrete symmetry

(19) |

where , etc. Note that Eq. (19) is not a gauge transformation, because does does not satisfy TBC [11]. This symmetry implies that the Polyakov lines in the twisted directions have zero expectation value, even in the perturbative phase of the theory.

We have done SU(3) simulations with twisted boundary conditions across two “spatial” boundaries and

(20) |

where , and across all three spatial boundaries , and

(21) |

We refer to these two cases as T and T boundary conditions, respectively.

Explicit representations of the twist matrices are not needed since, as shown in Ref. [11], one can absorb them by a field redefinition of the link variables, leaving only phase factors and multiplying the plaquettes at the corners of the twisted planes. We have done simulations with this field redefinition, and also using an explicit representation of the boundary conditions Eq. (14) using the matrices (for SU(3))

(22) |

and

(23) |

The boundary conditions Eq. (14) lead to an unusual quantization of the lattice momentum modes, as well as removing the zero modes. Making the usual substitution

(24) |

the boundary conditions take the form . A twisted plane wave basis is used to Fourier analyze the fields

(25) |

where, in order to obey the boundary conditions, the matrices must satisfy the algebra

(26) |

and where enforces a constraint on the mode sum, to be developed below (see Eqs. (30) and (31)).

The quantization conditions follow by iterating Eq. (26) times and using Eq. (16). One finds that momenta in twisted directions are quantized as if the SU() fields live on a lattice of length , rather than the actual length (although some modes are excluded)

(27) |

The extra momentum degrees-of-freedom come about because the color structure of each mode is unique, up to a phase. Substituting into Eq. (26) one finds, with a convenient choice of phase [11],

(28) |

These matrices are orthonormal under the trace

(29) |

Since the fields must be traceless one finds that a set of modes, including the zero modes, are excluded

(30) |

In the case of T boundary conditions, a further constraint emerges from Eq. (26)

(31) |

Hence there is a factor of more momentum modes with TBC, each of which has a single color degree-of-freedom, which is exactly the number of independent colors that one has for each momentum mode with PBC.

At tree-level the twisted gluon propagator in momentum space has the structure [11]

(32) |

where

(33) |

For the Wilson action in Lorentz gauges one has

(34) |

with and .

### iii.2 Suppression of phases

The Polyakov line along an untwisted direction is an order parameter for the degenerate vacua of the lattice theory, which correspond to the invariance of the Wilson action under the transformation

(35) |

The Polyakov line is also sensitive to the formation of domains between different phases. These nonperturbative effects must be suppressed if one is to use Monte Carlo simulations at large to extract perturbative quantities. In particular we will use simulation results for the Polyakov line itself to obtain the perturbative self-energy of a static quark.

In order to suppress nonperturbative phases we start the simulation with all links initialized to a “cold start,” . The probability of making a transition to an nontrivial phase in a local updating algorithm can then be reduced by working at sufficiently large lattice volumes. In fact our results for Wilson loops on lattices with PBC, presented in Sect. II, are in excellent agreement with finite-volume perturbation theory. However we find that nonperturbative phases are generated frequently on small lattices when PBC are used, and this occurs even at extremely large . On the other hand, we find that using TBC leads to a dramatic suppression of these effects compared to PBC, on lattices of the same size.

We illustrate the effects of phases with simulation results for lattices at . We show run time histories and scatter plots of the real and imaginary parts of the Polyakov line along an untwisted direction (hereafter taken to be the “temporal” direction ), where

(36) |

for a lattice of volume . Results for PBC are shown in Fig. 5, for T boundary conditions in Fig. 6, and for T boundary conditions in Fig. 7. We see that nonperturbative phases and domains render simulations with PBC useless for extracting perturbative quantities on small lattices. We also see that twisted boundary conditions create a barrier between phases, and that transitions between these phases and are essentially eliminated with T boundary conditions (with no tunneling events observed in millions of updates in the range of values considered here). In Sect. IV we show that the remaining finite-volume effects on lattices with T boundary conditions are very well described by perturbation theory for , even on volumes as small as .

## Iv Static-quark self-energy

### iv.1 Perturbation theory

In this section we consider the perturbative expansion of the self-energy of a static quark. We extract the self-energy from the gauge-invariant Polyakov line along an untwisted direction, which describes the propagation of a static quark across the entire time-extent of the lattice. One could also obtain the self-energy from large Wilson loops [15, 19], or from the gauge-fixed static-quark propagator [5]. This study however represents a good prototype for calculations of other more realistic perturbative quantities, such as quark masses.

We first define the self-energy on a finite lattice according to

(37) |

where, for comparison with our simulation results in the next subsection, we consider lattices with equal lengths along all sides. One then obtains the infinite-volume self-energy by taking the limit,

(38) |

We analyze the tadpole-improved self-energy. This is obtained by dividing the links in the Polyakov line by a mean field , . Hence the tadpole-improved self-energy is related to the unimproved self-energy by the addition of . We use the average plaquette to estimate the mean-field:

(39) |

The expansion of the self-energy to second order was computed in perturbation theory according to Eqs. (37) and (38) by Heller and Karsch, for an expansion in the bare coupling, with the result [14, 15]

(40) |

for the unimproved self-energy. Hereafter we consider only the tadpole-improved self-energy, which we denote by . We convert Eq. (40) to an expansion in the renormalized coupling at the appropriate scale using Eq. (8):

(41) |

where Eq. (5), with couplings evolved to , provides the tadpole subtraction.

In the next subsection we will compare Monte Carlo data for the self-energy with results from analytical perturbation theory on finite lattices. In order to extract the infinite-volume self-energy from the Monte Carlo simulations we must also make an extrapolation of measurements of done on finite volumes. We can gain some insight into the nature of the perturbative finite-volume corrections from some analytical considerations.

We define perturbative coefficients on a finite lattice according to

(42) |

For TBC the first-order term is given by

(43) |

to be compared with a calculation for PBC, ignoring the contribution from zero modes

(44) |

The constant that is subtracted from the momentum sums in the above expressions is the value of the one-loop tadpole-improvement counterterm, neglecting its very weak dependence on the lattice volume. We remind the reader that the mode sums in Eqs. (43) and (44) are different, due to the different quantization of the momentum components along the twisted and periodic directions. The two sets of boundary conditions yield identical results in the infinite-volume limit, where the color factor emerges in the case of TBC because averages to over infinitesimal momentum intervals.

Results for for the three boundary conditions are presented in Fig. 8, which shows that finite-volume effects are reduced with TBC, as suggested by Eq. (27). As is evident from the plots, finite-volume corrections are very well parameterized by a simple linear form in ,

(45) |

where . One can evaluate numerically from Eqs. (43) and (44) with the results (PBC), (T), and (T).

We expect the finite-volume correction to run with a coupling , evaluated at an infrared scale that is set by the box size:

(46) |

A physical interpretation of this functional form for is that the static quark experiences a perturbative Coulomb interaction with its images in the walls of the lattice. The different values of the coefficient for different boundary conditions also has a natural interpretation in this picture: TBC reduce finite-volume effects by effectively putting the image charges further away from the source charge.

Having established Eq. (46), one can deduce logarithms in in the self-energy at higher orders. For example, at second order one has

(47) |

where , and

(48) |

This follows from an expansion of the running coupling in Eq. (46), to second order in the coupling at a reference scale, such as . One can explicitly isolate the logarithm in the second order coefficient using existing perturbative calculations, which were done long ago by Heller and Karsch in the case of PBC [14], and which have also recently been done in Ref. [16] for TBC. Results of the perturbative calculations for the three boundary conditions are plotted in Fig. 9 over a range of lattice sizes. The dashed lines in Fig. 9 show fits to Eq. (47), where is constrained to the correct value; the fits are in excellent agreement with Eq. (48). Note that the curvature in the results for reveals the presence of the logarithm, particulary in the case of PBC and T boundary conditions, where the logarithm makes a significant contribution at the lattice sizes shown in Fig. 9.

In the next subsection we will use Eq. (46) to deduce the form of the logarithms in in the third order self-energy, which will help to constrain the infinite-volume extrapolation of the Monte Carlo data. We note that one should similarly be able to determine the leading logarithms in the finite-volume corrections to other quantities, which should likewise prove useful in Monte Carlo determinations of their perturbative expansions.

### iv.2 Self-energy from Monte Carlo simulations

We measured the static-quark self-energy in simulations done with T twisted boundary conditions at nine couplings. The simulation parameters are given in Table 5. Simulations were run on nine volumes , inclusive, at each of the nine couplings, for a total of 81 lattices. The number of measurements made on each volume, at all couplings except , were as follows: measurements for the lattices with inclusive, measurements for , for , for , for , and for (ten times as many measurements were made on each volume at ). One hundred configurations were skipped between measurements at all couplings, except at , where ten configurations were skipped between measurements. The observables were computed by binning the measurements in bin sizes of one hundred, which resulted in negligible autocorrelations at all couplings. The static energy and its error were computed from the binned ensembles using a standard jackknife analysis.

60.0 | 0.966311(1) | 0.008 | 0.00843 | |

23.8 | 0.913831(4) | 0.020 | 0.02338 | |

19.0 | 0.891415(5) | 0.025 | 0.03058 | |

16.0 | 0.869332(6) | 0.030 | 0.03824 | |

13.6 | 0.845310(7) | 0.035 | 0.04731 | |

12.0 | 0.822493(8) | 0.040 | 0.05676 | |

10.6 | 0.797025(10) | 0.045 | 0.06844 | |

9.5 | 0.771866(11) | 0.050 | 0.08138 | |

9.0 | 0.756142(13) | 0.053 | 0.09032 |

To demonstrate the reliability of the Monte Carlo method, we first use the simulation results to estimate the first- and second-order perturbative coefficients. In Fig. 10 we plot the quantity

(49) |

versus , for all values of . The dashed lines show the results of least-squares fits to Eq. (42), minimizing (Eq. (11)) using and for the first five orders in the expansion. The quality of the fits is very good in most cases, with -values typically in excess of about 20%, although the lowest -value in the fits is 3%.

We show the Monte Carlo results for in Fig. 11, where they are compared with finite-volume perturbation theory for T boundary conditions, Eq. (43). The data agree with perturbation theory within errors of only a few parts in .

Monte Carlo results for are shown in Fig. 12. In this case the fits to Eq. (42) were done with constrained to its perturbative value. We see that the Monte Carlo simulations also reproduce the results of second order T perturbation theory [16] over the full range of lattice sizes, within errors that are as small as a few parts in at several volumes.

The third-order term in the self-energy is not known from conventional perturbation theory. We use our simulation results to estimate , by redoing fits to Eq. (42) with both and constrained to their perturbative values. In order to determine the value of at infinite volume, we must also account for the systematic error due to the extrapolation from finite lattices. The leading finite-volume corrections come from logarithms in , which can be determined using renormalization-group arguments (cf. Eqs. (46)–(48)). Making use of these constraints considerably improves the accuracy of the extrapolation to infinite volume.

We show the Monte Carlo data for the third-order coefficient as a function of lattice size in Fig. 13, after subtracting the logarithms at ; the data are presented in terms of the residual

(50) |

The scale length in the leading logarithm is determined from second-order perturbation theory where, according to Eq. (47), . We evaluate the scale length from a fit to the T perturbation theory results illustrated in Fig. 9, which gives for the expansion in .

The Monte Carlo data are consistent with the logarithms in Eq. (50), which are found to dominate the extrapolation to the infinite-volume limit, even from these relatively small lattices. To extract the infinite-volume coefficient we fit the remaining finite-volume corrections to the form

(51) |

Figure 13 shows the results of a fit to Eq. (51) where is minimized using , and , for and (and for ). The fit yields

(52) |

Changing the order of the expansion in in Eq. (51) makes litle change in the fit value for . We conclude from these results that renormalized perturbation theory for the tadpole-improved self-energy is well behaved through third order, with the data in Fig. 13 clearly demonstrating that is of .

An estimate of the third-order term in the expansion of the unimproved self-energy, in the bare lattice coupling (Eq. (40)), has recently been reported using numerical simulations of the Langevin equations [19]. We can convert our result for from an expansion in to one in , using the inverse of Eq. (8). We find (without tadpole improvement), in agreement with the value reported in Ref. [19]. The bare coupling is clearly a very poor expansion parameter, with 96% of being absorbed by renormalization when a physical coupling is used.

We note that the simulations in Ref. [19] were done on much larger lattices than were used here. We were able to extract from smaller lattices because the leading finite-volume corrections were identified using renormalization-group methods. This determination of the third-order self-energy involved a modest computational effort. The entire set of simulations in the present analysis required only the equivalent of about 150 days of running on a single 1 GHz processor.

## V Summary and Outlook

The results presented here demonstrate that higher-order perturbative expansions are accessible in Monte Carlo simulations at large . An extensive theoretical analysis was presented together with the results of numerical simulations of a large set of Wilson loops and the static-quark self-energy. Twisted boundary conditions were investigated as a means of eliminating zero modes and suppressing nonperturbative finite-volume artifacts, and an extensive analysis of perturbative finite-volume corrections was made. Wilson loops provided a good quantity for a first test of the Monte Carlo method, since small loops are relatively insensitive to finite-volume effects. The calculation of the static-quark self-energy was considerably more involved, as it is very sensitive to nonperturbative phases and has large perturbative finite-volume corrections. The static-quark self-energy thus represents a good prototype for calculations of other more realistic perturbative quantities, such as quark masses.

The simulation results were shown to reproduce perturbation theory on finite lattices through second order to high precision, over a wide range of lattice sizes and couplings. Monte Carlo results for the fourteen smallest Wilson loops were found to agree with perturbation theory within the errors, with an accuracy of a few parts in for the first-order coefficients, and a few parts in for the second-order terms. The Monte Carlo results for the static-quark self-energy were found to agree with finite-volume perturbation theory over the full range of lattice sizes analyzed here, with an accuracy of a few parts in at first order, and a few parts in at second order. This precision was achieved with relatively little computational effort. New estimates of third-order terms for the Wilson loops and the static-quark self-energy were obtained to about 10% accuracy. Renormalization-group arguments were used to improve the quality of the extrapolation of the self-energy to infinite volume. The results demonstrate that renormalized perturbation theory for Wilson loops and the self-energy is well behaved through third order.

These methods can be directly applied to improved gluon actions, and can be extended to quark actions. We have done some work on large simulations for fermions in the nonrelativistic formulation of QCD, extending the preliminary studies reported in Refs. [1, 2]. We find that simulations of the additive energy and multiplicative mass renormalization reproduce results of one-loop perturbation theory, and can resolve the second-order terms in the expansion of these quantities, over a wide range of bare quark masses [5]. Further work in this direction is in progress.

###### Acknowledgements.

We are very grateful to Christine Davies for several crucial discussions. We thank Urs Heller for generously providing us with his programs for second-order coefficients. We also thank Richard Woloshyn, Mark Alford, Ron Horgan, Massimo DiPierro and Kent Hornbostel for helpful conversations. This work was supported in part by the National Science Foundation, the U.S. Department of Energy, and by the National Science and Engineering Research Council of Canada. HDT would also like to thank the United Kingdom Particle Physics and Astronomy Research Council, and the physics department of Cornell University, for support during part of this work.## Appendix:

Designing optimized Monte Carlo simulations

In this Appendix we consider the question of how to design the most efficient calculations through intelligent parameter choices, by using the techniques of constrained curve fitting [9]. As discussed in the Introduction, one objective of this analysis, in the context of short-distance Monte Carlo simulations like those in this paper, is to choose the couplings for the simulations so as to minimize the cost required to achieve a given precision in the perturbative coefficients.

As described in Ref. [9] and in Sect. II, the effects of truncation errors in fitting power series to Monte Carlo data may by estimated by augmenting with a function that constrains parameters, which are poorly determined statistically, to plausible values. Consider for example the logarithm of a Wilson loop, denoted by , which has the expansion

(53) |

We will use () to denote the set of couplings at which the simulations are done. We may define an augmented as in Eq. (11),

(54) |

where the second term on the right tends to constrain poorly determined parameters to the range , based on our prior experience with the power series. In the examples in this Appendix, we will use and .

A well designed calculation should minimize the errors in the final results, for a given amount of computer time. The uncertainties in the fit parameters are determined from the inverse of the Hessian matrix, which we denote by , where

(55) |

Then the uncertainty in is

(56) |

In the case of Eq. (53), an explicit expression for the Hessian matrix can be obtained

(57) |

where is the statistical error in . We note that a useful approximation for the statistical error in the Wilson loop (or its logarithm) is

(58) |

where and is independent of .

The optimal selection of the values at which the simulations are to be done may be determined by numerically minimizing the with respect to the . The optimal placement of ’s for the first couple of parameters may be guessed without doing much calculation. For example, since relative statistical errors are independent of , is obtained most accurately by running at the smallest possible coupling (avoiding round-off errors), thereby minimizing truncation errors from higher-order terms in the series expansion. The smallest total error in is then obtained by choosing a second coupling such that the statistical error in the slope of , which is , is equal to the truncation error, which is given by . Explicit numerical minimization reproduces expectations for these simple cases, but also produces optimized design parameters for arbitrary numbers of points and allocations of CPU time among them. The results of a three-point optimization are illustrated in Fig. 14. One sees that as the CPU time is increased, smaller and smaller couplings are obtained from the minimization calculations, though the optimal ’s do not fall quite as quickly as the fourth root of CPU time (as would be expected in a two-point fit).

Perturbative series in lattice QCD are typically used in nonperturbative Monte Carlo calculations with ’s in the range of about 0.15–0.30. To illustrate the effects of optimization calculations on the final result of a nonperturbative simulation, we consider the application of the perturbative series Eq. (53) to a simulation with . The results are shown in Fig. 15. The errors coming from the lowest-order coefficients are greatest when the short-distance simulations are done with low statistics, but these errors decrease most rapidly with CPU time; hence the order of the term that contributes the largest error to the nonperturbative quantity rises as a function of the CPU time of the perturbative simulations.

Similar results, different in detail, are obtained for other quantities. For example, for the third order coefficient for the Wilson loop, shown in Fig. 4, the statistical errors are and the truncation error is . The minimization formula gives the optimal placement for a single-point simulation as

(59) |

as expected.

Constrained curve fitting formulas for parameter fitting provide a concrete way of translating estimates of truncation errors into designs of efficient computations. For one or two parameters, they lead to exactly the same the same choices of parameters that intuitive guesswork provides. However, they also provide clear optimizations in the much more common situations in which we have too many parameters to guess about, or in which we are trying to design new runs to improve the results from imperfectly designed initial runs.

## References

- [1] W. B. Dimm, G. P. Lepage and P. B. Mackenzie, Nucl. Phys. B (Proc. Suppl.) 42, 403 (1995).
- [2] W. B. Dimm, Ph.D. thesis, Cornell University (1995).
- [3] W. B. Dimm, G. Hockney, G. P. Lepage and P. B. Mackenzie, unpublished [reported in C. T. H. Davies et al., Phys. Lett. B 345, 42 (1995); Phys. Rev. D 56, 2755 (1997).].
- [4] H. D. Trottier and G. P. Lepage, Nucl. Phys. B (Proc. Suppl.) 63, 865 (1998).
- [5] N. Shakespeare, Ph.D. thesis, Simon Fraser University (2000).
- [6] J. Juge, Nucl. Phys. (Proc. Suppl.) 94, 584 (2001).
- [7] G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 48, 2250 (1993).
- [8] C. T. H. Davies et al., Ref. [3].
- [9] G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel, P. B. Mackenzie, C. Morningstar, and H. D. Trottier, hep-lat/0110175.
- [10] G. ’tHooft, Nucl. Phys. B 153, 141 (1979).
- [11] M. Lüscher and P. Weisz, Nucl. Phys. B 266, 309 (1986).
- [12] A. Gonzalez Arroyo and C. P. Korthals Altes, Nucl. Phys. B 311, 433 (1988).
- [13] Some of this work was presented in preliminary form in G. P. Lepage, P. B. Mackenzie, N. H. Shakespeare and H. D. Trottier, Nucl. Phys. (Proc. Suppl.) 83, 866 (2000); and in Ref. [5].
- [14] U. M. Heller and F. Karsch, Nucl. Phys. B251 (1985) 254.
- [15] See also G. Martinelli and C. Sachrajda, Nucl. Phys. B 559, 429 (1999).
- [16] H. D. Trottier, G. P. Lepage, Q. Mason, M. A. Nobes, and K. Foley, in preparation. See also M. A. Nobes, H. D. Trottier, G. P. Lepage and Q. Mason, hep-lat/0110051.
- [17] See e.g. V. Gimenez, L. Giusti, G. Martinelli and F. Rapuano, JHEP 0003, 018 (2000); S. Collins, hep-lat/0009040; and Ref. [15].
- [18] A. Coste et al., Nucl. Phys. B 262, 67 (1985).
- [19] G. Burgio, F. Di Renzo, M. Pepe, and L. Scorzato Nucl. Phys. (Proc. Suppl.) 83, 935 (2000); F. Di Renzo and L. Scorzato, JHEP 0102, 020 (2001).
- [20] B. Allès et al., Phys. Lett. B 324, 433 (1994).
- [21] See e.g. I. Hinchliffe in Particle Data Group, D. E. Groom et al., Eur. Phys. J. C 15, 1 (2000).
- [22] M. Lüscher and P. Weisz, Phys. Lett. B 349, 165 (1995); Nucl. Phys. B 452, 234 (1995).
- [23] For this purpose we use a solution to the evolution equation that explicitly gives the dependence of the running coupling on the coupling at the reference scale. See, for example, G. Rodrigo, A. Pich, and A. Santamaria, Phys. Lett. B 424, 367 (1998); and G. Rodrigo and A. Santamaria, Phys. Lett. B 313, 441 (1993).
- [24] F. Di Renzo, E. Onofri and G. Marchesini, Nucl. Phys. B 457, 202 (1995).