Implementation of Symanzik's Improvement Program for Simulations of Dynamical Wilson Fermions in Lattice QCD

TL;DR

This paper details implementing Symanzik's improvement in lattice QCD simulations with dynamical Wilson fermions, showing minimal overhead and maintaining standard preconditioning techniques.

Contribution

It introduces a practical implementation of Sheikholeslami-Wohlert term within the Symanzik improvement framework for dynamical Wilson fermions in lattice QCD.

Findings

Overhead increases by about 20% when CG iterations exceed 200

Standard even-odd preconditioning remains effective

Implementation validated through analytical and numerical tests

Abstract

We discuss the implementation of a Sheikholeslami-Wohlert term for simulations of lattice QCD with dynamical Wilson fermions as required by Symanzik's improvement program. We show that for the Hybrid Monte Carlo or Kramers equation algorithm standard even-odd preconditioning can be maintained. We design tests of the implementation using analytically and numerically computed cumulant expansions. We find that, for situations where the average number of Conjugate Gradient iterations exceeds 200, the overhead is only about 20%.

Equations

$$\mathcal{Z}=\int \mathcal{D}U\mathcal{D}\bar{\psi }\mathcal{D}\psi \exp \left(-S_g-S_w-S_{sw}\right).$$

$$S_g=-\frac{\beta }{6}\sum _{P}Tr(U_P+U_P^{†}).$$

$$S_w=\sum _x\bar{\psi }(x)(D+m)\psi (x),$$

$$D$$

$${\nabla }_{\mu }\psi (x)$$

$${\nabla }_{\mu }^{\ast }\psi (x)$$

$$S_{sw}=\frac{i}{4}{c}_{sw}\sum _{x,\mu ,\nu }\bar{\psi }(x){\sigma }_{\mu \nu }{\mathcal{F}}_{\mu \nu }(x)\psi (x),$$

$${\mathcal{F}}_{\mu \nu }(x)$$

$$\mathcal{Z}$$

$$S_{eff}$$

$$Q(U{)}_{xy}$$

$$Q\equiv c_{0}\gamma_{5}\left(\begin{array}[]{cc}1+T_{ee}&M_{eo}\\ M_{oe}&1+T_{oo}\\ \end{array}\right)\;\;,$$

$$(T{)}_{xa\alpha ,yb\beta }=\frac{i}{2}{c}_{sw}\kappa {\sigma }_{\mu \nu }^{\alpha \beta }{\mathcal{F}}_{\mu \nu }^{ab}(x){\delta }_{xy}.$$

$$det(Q)$$

$$\hat{Q}$$

$$\mathcal{H}=\sum _{x,\mu }\frac{1}{2}Tr(H_{\mu }^2(x))+S_{eff}^{(eo)}(U_{\mu },{\varphi }^{†},\varphi ),$$

$$S_{eff}^{(eo)}$$

$$S_g[U_{\mu }]$$

$$S_{det}[U_{\mu }]$$

$$S_b[U_{\mu },{\varphi }^{†},\varphi ]$$

$$\varphi =\hat{Q}R,$$

$$U_{\mu }^{{}^{\prime }}(x)$$

$$i\delta H_{\mu }(x)$$

$$\delta S_{eff}^{(eo)}=\sum _{x\mu }Tr(F_{\mu }(x)\delta U_{\mu }(x)+F_{\mu }^{†}(x)\delta U_{\mu }^{†}(x)).$$

$$\delta S_b=-(X_o^{†}\delta \hat{Q}{Y}_o+Y_o^{†}\delta \hat{Q}{X}_o),$$

$$\delta \hat{Q}$$

$$\delta S_b=-({\hat{c}}_0/c_0)(X^{†}\delta QY+Y^{†}\delta QX),$$

$$X=\left(\begin{array}[]{c}-(1+T_{ee})^{-1}M_{eo}X_{o}\\ X_{o}\\ \end{array}\right)\;\;,$$

$$F_{\mu }^{(sw)}(x)=-\frac{{\hat{c}}_0{c}_{sw}\kappa }{8}(\text{mbox}\text{mainfont}diagramsinFig.2).$$

$$(■{)}_x=T{r}_{Dirac}(i{\gamma }_5{\sigma }_{\mu \nu }Y(x)\otimes X^{†}(x)+i{\gamma }_5{\sigma }_{\mu \nu }X(x)\otimes Y^{†}(x)).$$

$$\delta S_{det}=-2Tr(1+T_{ee}{)}^{-1}\delta T_{ee}.$$

$$(▲{)}_x=T{r}_{Dirac}[i{\sigma }_{\mu \nu }(1+T_{ee}(x){)}^{-1}].$$

$$F_{\mu }^{(det)}(x)=-\frac{\kappa c_{sw}}{4}(\text{mbox}\text{mainfont}diagramsinFig.3).$$

$$F_{\mu }(x)=V_{\mu }(x)+F_{\mu }^{(w)}(x)+F_{\mu }^{(det)}(x)+F_{\mu }^{(sw)}(x).$$

Full text

1 Introduction

In Wilson’s formulation of lattice QCD, when working at non-vanishing quark masses, chiral symmetry is violated. Its restoration requires the approach to the continuum limit and therefore small values of the lattice spacing $a$. Unfortunately, for small values of $a$ and light quark masses one needs large lattices and simulations easily become prohibitively costly. If, on the other hand, one stays at too large values of the lattice spacing, one has to face cut-off effects which are linear in the lattice spacing in the case of Wilson fermions. Indeed, in quenched simulations, using Wilson fermions, the lattice spacing effects were shown to be a severe problem [1].

A well-known remedy of the $O(a)$ effects is anchored in Symanzik’s improvement program [2] which leads to adding $O(a)$ correction terms to the lattice action. Such a term has been worked out for lattice QCD some time ago [3] and comes under the name of a Sheikholeslami-Wohlert (SW) term. It has a coefficient, $c_{sw}$, which has to be fixed so as to cancel the $O(a)$ effects. This kind of improvement has been used already in quenched simulations (see [4] for a review). Recently, using PCAC relations within the Schrödinger functional setup [1, 5], a non-perturbative determination of the coefficient was proposed and, at least in the quenched approximation, a significant effect of improvement could be verified.

It is a very natural next step to go beyond the quenched approximation and study the effects of improvement also in dynamical fermion simulations. In this paper we want to make a first step in this direction by describing how the SW-term can be implemented for simulations of lattice QCD using molecular dynamics algorithms like Hybrid Monte Carlo [6] or Kramers equation [7, 8]. An important observation is that for the implementation conventional even-odd preconditioning [9] can be maintained 111Of course, also different preconditioning techniques like the one proposed in [10] could be used..

We will give analytically and numerically computed cumulant expansions that can serve as checks of the implementation. Finally we will discuss the crucial question, whether the gain in reducing the lattice spacing effects will merit the overhead of adding a SW-term to Wilson’s fermion action.

The theory is established on a Euclidean space-time lattice with size $L^{3}\times T$. With lattice spacing set to unity from now on, the points on the lattice have integer coordinates $(t,x_{1},x_{2},x_{3})$ which are in the range $0\leq t<T;0\leq x_{i}<L$. A gauge field $U_{\mu}(x)\in SU(3)$ is assigned to the link pointing from point $x$ to point $(x+\mu)$, where $\mu=0,1,2,3$ designates the 4 forward directions in space-time. The gauge fields assume periodic boundary conditions. The quark fields $\psi_{Aa\alpha}(x)$ are defined on each site of the lattice, where $A$,$a$ and $\alpha$ are flavor, color and Dirac indices respectively. The partition function for Wilson QCD is given by,

[TABLE]

The first term $S_{g}$ is the pure gauge action and is given by

[TABLE]

The symbol $U_{P}$ represents the usual plaquette value on the lattice. The bare gauge coupling $g_{0}$ is related to $\beta$ by $g^{2}_{0}\equiv 6/\beta$. The second term in the action is the usual Wilson fermion action:

[TABLE]

where $m$ is the bare quark mass. The Wilson difference operator $D$ which appears in the above expression is given by:

[TABLE]

The third term $S_{sw}$ in the action is the Sheikholeslami-Wohlert (or clover) term. We write it as:

[TABLE]

where $c_{sw}$ is a parameter that has to be tuned such that this term will eliminate the $O(a)$ effects in the on-shell physical observables (Symanzik improvement). To tree level, this parameter should be set to unity. The Dirac matrices $\sigma_{\mu\nu}$ are defined in the usual way via $\gamma$-matrices. The antisymmetric and antihermitian tensor ${\cal F}_{\mu\nu}(x)$ is a function of the gauge links and given by:

[TABLE]

We can represent this quantity by its corresponding diagram (see Fig. 1), which suggests the name “clover”.

We will assume that the mass matrix in equation (3) is proportional to the unit matrix in flavor space, namely we are considering two flavors of Wilson fermions with degenerate masses. The fermion fields in the partition function can be integrated out, resulting in the fermion determinant. This determinant is then “integrated in”, using the pseudofermion (boson) fields $\phi^{\dagger}(x)$ and $\phi(x)$ for the convenience of the simulation. We then write the partition function and effective action of the theory as:

[TABLE]

where the matrix $Q$ is related to the fermion matrix $M$ by $Q=c_{0}\gamma_{5}M$ and, including the SW-term, is given by:

[TABLE]

with $\kappa=(8+2m)^{-1}$ and $c_{0}=(1+8\kappa)^{-1}$. The notation $\sigma_{\mu\nu}{\cal F}_{\mu\nu}$ in the above expression, and henceforth, implies the summation over $\mu$ and $\nu$.

2 Implementation of the SW-term

In this section we will discuss how the SW-term is implemented for simulations using Hybrid Monte Carlo or Kramers equation algorithms. Although these two algorithms show a comparable performance, the Kramers equation algorithm seems to be advantageous due to possible problems of non-reversibility in the Hybrid Monte Carlo algorithm [8]. Since the use of preconditioning techniques in fermion simulations appears to be very important we will discuss directly how to install the SW-term with, in our case, even-odd preconditioning, see also [11] for a similar discussion. Let us write the matrix $Q$ in equation (8) as

[TABLE]

where we have introduced the matrix $T_{ee}$($T_{oo}$) on the even (odd) sites as

[TABLE]

The off-diagonal parts $M_{eo}$ and $M_{oe}$, which connect the even with odd and odd with even lattice sites respectively, are just the conventional Wilson hopping matrices. Preconditioning is now realized by writing the determinant of $Q$, apart from an irrelevant constant factor, as

[TABLE]

The constant factor $\hat{c}_{0}$ is given by $\hat{c}_{0}=1/(1+64\kappa^{2})$. For the simulation with even-odd preconditioning, one introduces the Hamiltonian:

[TABLE]

where $H_{\mu}(x)$ represents the momentum conjugate to the gauge field $U_{\mu}(x)$ and takes values in the Lie algebra of $SU(3)$, i.e. a traceless hermitian $3\times 3$ matrix. The effective action $S^{(eo)}_{eff}$ for the even-odd preconditioned case is given by:

[TABLE]

The first term $S_{g}$ in $S_{eff}^{(eo)}$ is just the usual Wilson pure gauge plaquette action. The second term involves the determinant of the matrix $(1+T_{ee})$ on the even lattice sites. We shall refer to it as the determinant contribution. The third term is the determinant of the preconditioned matrix $\hat{Q}$ which acts only on the boson fields at odd lattice sites. We will call $S_{b}$ the bosonic contribution.

The preconditioned matrix $\hat{Q}$ contains the inverse matrix $(1+T_{ee})^{-1}$ which is needed only locally at a given lattice point $x$. Separating the lower and upper components of the fermion fields, the computation of the inverse amounts to the inversion of two complex $6\times 6$ matrices per lattice site. This can be done using, for example, a Householder triangularization [12] which is advantageous on parallel computers, especially on the Alenia Quadrics (APE) machines we are using here [13]. In order to study, whether there might occur problems with inverting the matrix $(1+T_{ee})$, we measured the lowest value $D_{min}$ of the determinants $det(1+T_{ee}(x))$ for every configuration. As an example, for $\beta=5$, $\kappa=0.150$ and $c_{sw}=1.52$ we find that $D_{min}=0.584(2)$. Noting that the critical hopping parameter $\kappa_{c}$ is lowered for SW-improved fermions as compared to conventional Wilson fermions [14] and, since for larger values of $\beta$ and(or) lower values of $\kappa$, $D_{min}$ is increasing, we do not expect problems with the matrix inversion for realistic QCD simulations.

Since the bosonic part $S_{b}$ is quadratic in the $\phi$ fields, they are generated at the beginning of each molecular dynamics trajectory via

[TABLE]

where $R$ is a random spinor field taken from a Gaussian distribution of norm one. The gauge fields and its corresponding momenta are updated according to:

[TABLE]

where $[\cdots]_{T.A.}$ stands for the traceless antihermitian part of the matrix [15] and the quantity $[U_{\mu}(x)F_{\mu}(x)]_{T.A.}$ is the “total force” associated with the link $U_{\mu}(x)$ which will be discussed in more detail below. In the simulation the update eq.(2) is realized via a leapfrog integration scheme.

In molecular dynamics algorithms, $F_{\mu}$ is the “coefficient” in the change of the effective action when an infinitesimal change of the gauge link $\delta U_{\mu}(x)$ is applied, i.e.

[TABLE]

The quantity $F_{\mu}(x)$ receives three contributions. The first one, $V_{\mu}(x)$, is coming from the pure gauge part. The second is from the determinant of the SW-term on the even sites. The third one is coming from the change in $\phi^{\dagger}\hat{Q}^{-2}\phi$. We will now discuss the last two terms in some detail.

Let us start with the third contribution, namely the variation of the effective action from the bosonic term $S_{b}$. A direct variation of this term gives:

[TABLE]

where $X_{o}\equiv\hat{Q}^{-2}\phi_{o}$ and $Y_{o}\equiv\hat{Q}^{-1}\phi_{o}$ are fields defined only on odd lattice sites. They are obtained by inverting the preconditioned matrix $\hat{Q}^{2}$ with a given source $\phi_{o}$.

Using

[TABLE]

we find that the variation $\delta S_{b}$ can be written as

[TABLE]

where now, the fields $X$ and $Y$ are defined over the full lattice. For the field $X$, we have:

[TABLE]

and similarly for the field $Y$. The variation of the original matrix $Q$ again consists of two parts. The contribution from the off-diagonal elements is the same as in the conventional Wilson case. This term will be denoted as $F^{(w)}_{\mu}(x)$ for a given link $U_{\mu}(x)$. Therefore, the only new ingredient is from the diagonal terms which explicitly involves the contributions of the SW-term. To obtain this contribution, let us focus on the variation of one particular link variable $U_{\mu}(x)$, keeping all others constant, and try to collect all possible contributions. One easily finds out that, in order for the “clover leaves” of the SW-term to contain the chosen link $U_{\mu}(x)$, the center of the clover can only be at lattice points $x$, $x+\mu$, $x\pm\nu$ and $x+\mu\pm\nu$, where $\nu$ is some direction different from $\mu$.

After realizing this, it is just a matter of book-keeping to gather all the possible terms. We will represent them collectively in a diagrammatic fashion, see Fig. 2. We write the contribution, $F^{(sw)}_{\mu}(x)$, from these diagrams as:

[TABLE]

The diagrammatic representations have the following simple interpretation. Starting from point $x+\mu$, one can, following the direction of the arrows on the lines, go around the plaquette and finally arrive at point $x$. In doing this, one performs matrix multiplications along the way. Each line segment with an arrow on it represents the corresponding gauge link at that point. A black square at a given point, with a particular choice of $\mu$ and $\nu$, stands for an “insertion” of a $3\times 3$ matrix in color space. This matrix is given by the following expression:

[TABLE]

The trace in the above formula, as indicated, is done only in Dirac space and $X\otimes Y^{\dagger}$ stands for a direct product in color space. For a fixed $\mu\nu$ pair, one can take two different paths to arrive at point $x$, namely the “upward” path, which passes through point $(x+\nu)$, and the “downward” path, which passes through point $(x-\nu)$. The path that winds downward gets a relative minus sign for its contribution. Finally, one should sum up all the contributions from all possible $\nu\neq\mu$ values.

Let us turn now to the determinant contribution from the even sites. We have, when taking the variation of the action,

[TABLE]

This contribution can also be easily represented by a set of diagrams. In fact, by almost the same set of diagrams as

in the previous case, except for the following changes. First the black square insertions at various points are replaced by the black triangles at the same points. And each black triangle, at some chosen point $x$ and direction $\nu\neq\mu$, represents the following expression:

[TABLE]

The second modification is that these insertions should only be done at even lattice points, otherwise the contribution is set to zero. With these modifications, the contribution from the determinant $F^{(det)}_{\mu}(x)$ is given by

[TABLE]

Finally, $F_{\mu}(x)$ is written as a sum of all the contributions as:

[TABLE]

Let us mention, that even-odd preconditioning can also be achieved in a manner slightly different from the case just discussed. In this case one treats the even and odd points more symmetrically. This amounts to rewriting the determinant of the fermion matrix as:

[TABLE]

Then, the determinant of $\hat{M}^{\dagger}\hat{M}$ is realized by introducing the bosonic Gaussian fields while the determinants of $(1+T_{ee})$ and $(1+T_{ee})$ could be treated the same way as discussed before. In this version the determinant contributions are then summed up on all lattice sites, even and odd. At present, it is not clear, whether one of the two possibilities of preconditioning the fermion matrix should be preferred.

3 Tests of the implementation

Applying an infinitesimal change to the gauge links, the subroutine which evaluates $F_{\mu}(x)$ was first checked using eq.(16) by computing the change of the action explicitly and by using $F_{\mu}(x)$. Besides this, the focus of our additional tests will concern the SW-term alone. For this purpose, it is advantageous to study cases in which the effects of the SW-term have been singled out. The SW-term that enters the simulation program only depends on the parameter combination $\bar{c}_{sw}\equiv c_{sw}\kappa$. Therefore, one can take the limit $\kappa\rightarrow 0$ but keep $\bar{c}_{sw}$ at some prescribed value. In this way, the effects of the Wilson hopping matrices are completely switched off and the fermionic part of the action depends only on the SW-term. Of course, the pure gauge parameter $\beta$ is still at our disposal and the SW-term can be tested at various values of $\beta$.

3.1 Tests at $\beta=0$

As a first test, we study the model in the limit of $\beta=0$ and $\kappa=0$, however, $c_{sw}\kappa$ is kept finite. In this limit, the effective action of the theory is simply:

[TABLE]

with $T(x)$ as given in eq.(10). The advantage of taking this limit is that, in this simplified case, the average plaquette value can be evaluated analytically as a power series expansion in terms of $\bar{c}_{sw}$. For the expectation value of a plaquette, $P\equiv(1/6)tr(U_{P}+U^{\dagger}_{P})$, to lowest order one gets

[TABLE]

To work out the next leading correction of this expansion requires

some more work. Basically one is led to the following general expression:

[TABLE]

with $<O>_{0}=\int dUO(U)$ and $dU$ being the normalized group measure. In the above formula, we have used the shorthanded notation $\sigma{\cal F}$ to stand for $\sigma_{\mu\nu}{\cal F}_{\mu\nu}/2$. The last two terms correspond to the connected contribution of the operators $P$ and $Tr(\sigma{\cal F})^{2}$. The calculation of these terms is quite straightforward but rather technical. We list some of the details in the appendix. This calculation results in:

[TABLE]

We made simulations with SW-improved action in the Schrödinger functional setup 222For the timelike plaquettes touching the boundary within the Schrödinger functional, eq.(31) is replaced by $-(1/24)\bar{c}_{sw}^{2}(1+(491/384)\bar{c}^{2}_{sw})$. at various values of $\bar{c}_{sw}$ at $\beta=0$ on a $4^{4}$ lattice and the results are shown in Fig. 4 together with the prediction of eq. (29) and eq. (31). For small enough values of $\bar{c}^{2}_{sw}$ (less than $0.01$), the lowest order of the expansion is already satisfactory. For larger $\bar{c}^{2}_{sw}$ values, the deviation of the data from the lowest order prediction is quite significant. However, if one also includes the second order contribution, all data points are in good agreement with the anticipated values up to $\bar{c}^{2}_{sw}\sim 0.1$.

3.2 Tests at nonvanishing $\beta$

For a nonvanishing value of $\beta$ and any observable $O$, one can obtain an expansion in terms of numerically computed cumulants at $\bar{c}_{sw}=0$. We have

[TABLE]

In the above formula, the notation $<...>_{0}$ stands now for the Monte Carlo average at $\bar{c}_{sw}=0$. We have performed a high statistics run at a fixed $\beta$ value with $\bar{c}_{sw}=0$, measured the cumulants above and constructed the expansion. On the other hand, we have made simulations at a nonvanishing value of $\bar{c}_{sw}$ directly and compared with the cumulant expansion. This test was performed at $\beta=6$ on a $4^{4}$ lattice and the result is illustrated in Fig. 5.

Obviously, for the two quantities that have been measured, namely the average plaquette value and the average value of $Tr({\cal F}_{\mu\nu}{\cal F}_{\mu\nu})$, the measured values from the simulation agree very well with the prediction of the cumulant expansion. The width of the shaded bands in the figure represents the statistical uncertainty of the cumulant expansion which is constructed from a run at $\bar{c}_{sw}=0$ with $5.5$ million measurements. In this analysis, we only took into account the blocked statistical errors of the run at $\bar{c}_{sw}=0$. The $O(\bar{c}^{3}_{sw})$ term in the cumulant expansion (32) is also included, the effect of which is however rather small. It is clear from the figure that the systematic errors in this case should be smaller compared with the statistical errors for the parameter range we are considering here. Of course, by measuring more cumulants at $\bar{c}_{sw}=0$, one could also construct higher order corrections in the cumulant expansion.

4 Timing of the Program

When comparing with the conventional Wilson QCD simulation, the timing of the version with the improved fermion action is a crucial piece of information.

In a typical QCD simulation using molecular dynamics algorithms, the conjugate gradient (CG) iterations dominate the CPU time of the program. Therefore, in the limit of number of conjugate gradient iterations $N_{CG}$ going to infinity, the time of the matrix-vector multiplication is the most crucial part. In our version of Symanzik improvement, this multiplication turns out to be only slightly slower than the conventional Wilson case by about 15 percent. The force evaluation, in this case, adds some overhead to the program which does not depend on $N_{CG}$. It is convenient to express all times in units of the matrix-vector multiplication of the conventional Wilson case.

To be specific, we propose the following formula for the time of the program:

[TABLE]

In the above formula, the coefficient $c^{(2)}_{w}=2$ is, by definition, the number of matrix-vector multiplications needed for each CG iteration in the conventional Wilson case. The coefficient $c^{(1)}_{w}$ represents the cost of other operations for each CG iteration (linear combinations, inner products etc.). Typically, this coefficient is quite small. The coefficient $c^{(0)}_{w}$ is the overhead that does not depend on $N_{CG}$. The value of $c^{(0)}_{w}$ depends on the implementation of the program. However, in the asymptotic region where $N_{CG}$ is large, the effect of $c^{(0)}_{w}$ becomes irrelevant.

The quantities $c^{(0)}_{sw}$, $c^{(1)}_{sw}$ and $c^{(2)}_{sw}$ are the corresponding coefficients for the program with the Sheikholeslami-Wohlert action. In this case, $c^{(2)}_{sw}=2.3$ is slightly larger than the corresponding value in the Wilson case. The value of $c^{(1)}_{sw}=c^{(1)}_{w}$ remains the same and we find, for our implementation of the SW-term on the APE computer, the coefficient $c^{(0)}_{sw}\sim(50+c^{(0)}_{w})$ to be larger than that of the Wilson case. In Fig. 6 we plot the CPU time of our program with improved fermions relative to the conventional Wilson case as a function of the number of CG-iterations per molecular dynamics step. The solid and dashed curves in the figure correspond to the choice of $c^{(0)}_{w}=0$ and $c^{(0)}_{w}=50$ respectively, which we consider to be two rather extreme cases. For typical simulations, the curve should lie between these two. In any case, we see that, for the most interesting physical situation in which $N_{CG}\sim 200$ or more, the simulation with improved fermions is only about 20% slower, quite independent of the value of $c^{(0)}_{w}$.

5 Conclusions

We demonstrated in this paper that the SW-term as required by Symanzik’s improvement program to cancel $O(a)$ effects in Wilson QCD can be implemented straightforwardly for simulation algorithms like Hybrid Monte Carlo or the Kramers equation. For the implementation even-odd preconditioning can be maintained which is an important aspect concerning the performance of the algorithms. We designed tests for checking the implementation of the SW-term by analytical and numerical computations of cumulant expansions.

The most important outcome of our work is illustrated in Fig. 6 which compares the time of a program with conventional Wilson fermions with the one where the SW-term has been added. The figure nicely demonstrates that the overhead due to the SW-term in the lattice action is only about 20$\%$ when the number of conjugate gradient iterations exceeds about 200 iterations, which is a number easily reached in realistic simulations of lattice QCD. On the other hand, Symanzik improvement will enable us to perform QCD simulations with much smaller lattice artifacts than the conventional Wilson fermions. Since the cost of a QCD simulation is proportional to a large power of $a^{-1}$, any small factor of improvement in lattice spacing will reduce the cost of the simulation significantly.

This brings us to the conclusion that simulations with SW-improved dynamical Wilson fermions can be done without much more cost than conventional Wilson fermion simulations. We expect that the cancellation of the strong $O(a)$ effects, as have been seen in quenched simulations, will by far merit the small overhead due to the addition of the SW-term.

A last remark concerns the multiboson technique for simulating lattice QCD [16, 17]. This relatively new method shows a comparable performance to the molecular dynamics algorithms and can be considered as a real alternative to the Hybrid Monte Carlo or Kramers equation algorithms. However, for the multiboson technique it is not straightforward to implement the SW-term, if one wants to use heatbath or over-relaxation algorithms. We propose therefore, to use exact local versions of the Hybrid Monte Carlo or Kramers equation algorithms [18] for the simulation. For these algorithms, the implementation of the SW-term described here directly applies.

6 Acknowledgements

We would like to thank M. Lüscher and U. Wolff for very helpful discussions. We thank P. Weisz and U. Wolff for offering their matrix multiplication routines for the preconditioned fermion matrix to us. The numerical simulations have been performed on the APE/Quadrics computers at DESY-IFH (Zeuthen). We thank the staff of the computer center at Zeuthen for their support.

Appendix A Analytic computation of the cumulant expansion at $\beta=0$

In this appendix, we discuss the analytic calculation of the cumulant expansion for the average plaquette value at $\beta=0$ in some detail. Our aim is to calculate the ensemble average of the plaquette which is defined as $P=(1/6)tr(U_{P}+U^{\dagger}_{P})$. We will only consider the case of periodic boundary conditions for the gauge fields here. The plaquette is specified by its (geometrical) “orientation”, which designates the plane that the plaquette lies in, and its two possible “circulations”, one corresponding to $U_{P}$ and the other to $U^{\dagger}_{P}$. When we expand the effective action (28) into power series of $\bar{c}_{sw}$, at each order, we will have products of plaquettes of various orientations and circulations. However, in order to obtain a nonvanishing contribution after group integrations, all the plaquettes that appear in the product have to be “tiled” according to the following simple rule:

• For any plaquette appearing in the expansion, the number of factors of one circulation minus the number of factors of the opposite circulation has to be [math] (mod 3).

The average plaquette value, up to $O(\bar{c}^{4}_{sw})$, is given by equation (30) in the main text. To obtain the contribution of $O(\bar{c}^{4}_{sw})$, one has to evaluate the contributions from the following two cases: (1) the connected contributions of $P$, $Tr(\sigma{\cal F})^{2}$ and $Tr(\sigma{\cal F})^{2}$ and (2) the contributions from $PTr(\sigma{\cal F})^{4}$. We will now study these two cases separately. Since $tr(U_{P})$ and $tr(U^{\dagger}_{P})$ always give the same contribution, we will just consider the former. In the calculation below, the following common factor always appears:

[TABLE]

which comes from the definitions of the plaquette and ${\cal F}_{\mu\nu}$.

A.1 The connected contribution

$<PTr(\sigma{\cal F})^{2}Tr(\sigma{\cal F})^{2}>_{c0}$

The connected contribution $<PTr(\sigma{\cal F})^{2}Tr(\sigma{\cal F})^{2}>_{c0}$ is defined to be

[TABLE]

On the right hand side of the above equation, the first term receives contributions of the type $tr(U_{P})tr(U_{P}U_{P})tr(U_{P^{\prime}}U^{\dagger}_{P^{\prime}})$ which exactly cancels the second term. The remaining contributions only come in from two different cases which we will call (a) and (b). In case (a), the two factors of $Tr(\sigma{\cal F})^{2}$ are sitting at neighboring corners of the plaquette to be measured. The link which connects these two corners is also an edge of another plaquette $P^{\prime}\neq P$ which has the same geometrical orientation. Then, the combination of the type $tr(U_{P})tr(U_{P}U_{P^{\prime}})tr(U_{P}U^{\dagger}_{P^{\prime}})$ fulfills the rule and gives a contribution of

[TABLE]

In case (b), the situation is the same as in case (a) except that two factors of $Tr(\sigma{\cal F})^{2}$ are sitting at the same corner of the plaquette to be measured. Now, since we have $3$ different ways of choosing the other plaquette $P^{\prime}$, we get a contribution:

[TABLE]

This concludes the calculation of the contribution $<PTr(\sigma{\cal F})^{2}Tr(\sigma{\cal F})^{2}>_{c0}$.

A.2 The contribution from $<PTr(\sigma{\cal F})^{4}>_{0}$

In this case, there could be two possibilities concerning the Dirac structure of the factors from $\sigma{\cal F}$:

• (A) All $4$ factors have the same geometrical orientation as the plaquette to be measured.
• (B) Two of the $4$ factors have the same geometrical orientation as the plaquette to be measured, the other two take orthogonal orientations.

We will now discuss these contributions separately.

For case (A), we have the following $3$ possibilities which we will list as (A1), (A2) and (A3). In all these cases, the Dirac trace simply gives a factor of $4$. In case (A1), we consider the combination $tr(U_{P}U_{P}U_{P^{\prime}}U^{\dagger}_{P^{\prime}})$ , $P\neq P^{\prime}$. Note that there are $4$ such terms in $Tr(\sigma{\cal F})^{4}$. The contribution is:

[TABLE]

In case (A2), the situation is almost the same as in case (A1) with the exception that the structure is of the type $tr(U_{P}U_{P^{\prime}}U_{P}U^{\dagger}_{P^{\prime}})$. Since there are $2$ such terms in $Tr(\sigma{\cal F})^{4}$, the contribution is:

[TABLE]

In case (A3), at least one of the $\sigma{\cal F}$ factors contributes a $U^{\dagger}_{P}$ which has a circulation opposite to that of the plaquette to be measured. This contribution is given by:

[TABLE]

where the group integral $I$ is given by

[TABLE]

This integral, as well as the other group integrals appearing in the different contributions, can be evaluated using the techniques described in appendix B. The value of $I$ in (41) equals $5$. This concludes the discussion of case (A).

For the case (B), we have the following two possibilities which we call (B1) and (B2).

In case (B1), the situation is very similar to case (A1) above, i.e. the contribution has the form $tr(U_{P}U_{P}U_{P^{\prime}}U^{\dagger}_{P^{\prime}})$. In this case, the Dirac trace will just give a factor of $4$. The two factors $U_{P^{\prime}}$ and $U^{\dagger}_{P^{\prime}}$ can take the $5$ orientations orthogonal to the one specified by the measured plaquette. They can also take all $4$ possible leaves and $2$ possible circulations in the clover. This contribution turns out to be the dominant one at this order. We have:

[TABLE]

Case (B2) is closely related to (A2) and the structure is $tr(U_{P}U_{P^{\prime}}U_{P}U^{\dagger}_{P^{\prime}})$. However, since $P^{\prime}$ does not have the same geometrical orientation as $P$, the Dirac trace should be done with more care. Basically, for any given orientation of the measured plaquette, among the other $5$ orthogonal orientations, only one of them gives a Dirac trace of $4$ while the other $4$ orientations give a Dirac trace of $(-4)$ due to the anticommuting properties of the $\sigma$ matrices. Therefore, effectively there appear to be only $3$ orientations to contribute a Dirac trace $(-4)$. We get

[TABLE]

This concludes our calculation for case (B).

Adding all contributions up, gives the result quoted in the main text.

Appendix B Evaluation of group integrals for $SU(3)$

Here we discuss some techniques to evaluate the group integrals for the group $SU(3)$. The integrals that we wish to evaluate are of the form

[TABLE]

where $U$ is an element of $SU(3)$ and the measure $dU$ stands for the normalized Haar measure of the group. The first step is to use the Cayley relation for an $SU(3)$ matrix:

[TABLE]

which could be established easily from the fact that $(U-\lambda_{1})(U-\lambda_{2})(U-\lambda_{3})=0$, where $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are the eigenvalues of $U$ and unit matrices are to be substituted appropriately. With the help of this relation, we can transform terms with powers of $U$ higher than $1$ into powers of $1$, [math] and $(-1)$. This means that, without loss of generality, we can discuss the following integral:

[TABLE]

This integral can be done by realizing that the integrand above is nothing but the character of the representation $\mathbf{3}^{i_{1}}\otimes\bar{\mathbf{3}}^{j_{1}}$, assuming $U$ is in the $\mathbf{3}$ representation. Therefore, the above integral is equal to the number of trivial representations in this direct product. All one has to do now is to decompose the direct product $\mathbf{3}^{i_{1}}\otimes\bar{\mathbf{3}}^{j_{1}}$ and count how many times the trivial representation appears.

As an example, let us consider the integral

[TABLE]

Using Cayley relation, this is transformed into:

[TABLE]

The last term just gives $4$ and the third term gives $2$. One easily finds out that $\mathbf{3}\otimes\mathbf{3}\otimes\bar{\mathbf{3}}\otimes\bar{\mathbf{3}}$ contains $2$ trivial representations. So the second term yields $(-4)\cdot 2=(-8)$. The first term requires the decomposition of $\bar{\mathbf{3}}\otimes\mathbf{3}\otimes\mathbf{3}\otimes\mathbf{3}\otimes\mathbf{3}$ and one finds that it contributes a $3$. Therefore the integral finally yields a value of $1$.