text:
# 1 Introduction

## 1 Introduction

In this article we propose an algorithm for dynamical simulations of smeared link fermions. We consider a system described by the action

$$S=S_g(U)\mathrm{tr}\mathrm{ln}\left[Q^{}(V)Q(V)\right],$$
(1)
where $`S_g(U)`$ is the pure gauge action and $`Q(V)`$ is the fermionic matrix. The gauge connections between the fermions are smeared links $`V`$ which are constructed deterministically from the dynamical thin links $`U`$. Since each smeared link is a local combination of a finite number of thin links, the system where the fermions couple via smeared links is in the same universality class as the system with thin links. In the action we keep the smearing of the links unchanged as the continuum limit is approached.

Smearing of the gauge links in the action and operators of a theory is part of many improvement programs. It has been demonstrated that smeared links improve flavor symmetry with staggered fermions \[?,?,?\] and chiral symmetry with Wilson-type clover fermions \[?\]. They are useful in the construction of an overlap Dirac operator \[?,?,?\] and arise naturally in the framework of perfect actions \[?,?,?\]. As an example for improved operators we mention the classically perfect Polyakov loop in pure gauge theory \[?\]. In order to eliminate lattice artifacts the static source has to become an extended object. Wilson loops constructed from smeared links were found to improve the statistical accuracy of potential measurements in \[?,?\] also.

The problem with eq. (1) is that the standard simulation algorithms for dynamical fermions (Hybrid Monte Carlo (HMC) or R algorithms) involve the computation of the gauge force $`Q(V)/U`$, which is either very complicated or likely impossible if the smeared links are projected onto SU(3). Present large scale dynamical simulations with standard algorithms \[?,?\] use one level of non–projected smearing. Quenched studies indicate that it would be desirable to simulate dynamical fermions with smoother links, either obtained by a more complex and effective smearing procedure like the hypercubic (HYP) blocking \[?\] or by iterative smearing with projection onto SU(3) after each smearing step.

## 2 The HYP action

The definition of HYP links and their properties are discussed in \[?,?\]. In the natural formulation of four flavors of staggered quarks, the continuum SU(4) flavor symmetry is broken to U(1) at finite lattice spacing. The pion spectrum has only one true Goldstone pion $`\pi _\mathrm{G}`$, the other 14 pions remain massive when the quark mass approaches zero. The mass splitting between a non–Goldstone pion $`\pi `$ and $`\pi _\mathrm{G}`$ can be parametrized by the quantity \[?\]

$$\delta _2=\frac{m_\pi ^2m_\mathrm{G}^2}{m_\rho ^2m_\mathrm{G}^2}$$
(2)
evaluated at a fixed $`m_\mathrm{G}/m_\rho `$ ratio. With a flavor symmetric action $`\delta _2=0`$ for all the pions.

We computed $`\delta _2`$ at $`m_\mathrm{G}/m_\rho =0.55`$ for the two lightest non-Goldstone pions, $`\pi _{i,5}`$ with flavor structure $`\gamma _i\gamma _5`$ and $`\pi _{i,j}`$ with flavor structure $`\gamma _i\gamma _j`$, on a set of quenched $`8^3\times 24`$ lattices generated with Wilson pure gauge action at $`\beta =5.7`$ ($`a=0.17\mathrm{fm}`$), using three different valence quark actions: standard staggered action with thin or HYP smeared links and improved Asqtad action of the MILC collaboration \[?\]. The results are shown in table 1. Flavor symmetry violations with the HYP action are reduced by close to an order of magnitude with respect to the thin link action and by a factor of two with respect to the Asqtad action.

The HYP blocking was designed to improve flavor symmetry by reducing the fluctuations of the gauge links that couple to the different flavor and Dirac components of staggered fermions inside a hypercube. The parameters of the HYP transformation have been optimized non–perturbatively to minimize the short distance fluctuations of the smeared gauge field \[?\]. A first non–trivial check of the HYP construction comes from perturbation theory: The cancellation of the flavor symmetry violating terms at tree level in the HYP staggered action gives consistent results for the HYP parameters \[?\]. We looked further at the scaling properties of the HYP action. Figure 1 shows the rho mass as a function of the lattice spacing squared, all in units of $`r_0`$ \[?\]. In order to compare a result from a dynamical HYP simulation at gauge coupling $`\beta =5.1`$ and bare quark mass $`am=0.04`$, we chose to extrapolate the quenched data to $`m_\mathrm{G}/m_\rho =0.52`$, a slightly different value than the one used in table 1. The quenched results correspond to $`\beta =5.7`$ ($`r_0/a=2.96(2)`$) and $`\beta =6.0`$ ($`r_0/a=5.36(4)`$) with the Wilson gauge action and to $`\beta =8.0`$ ($`r_0/a=3.687(13)`$) and $`\beta =8.4`$ ($`r_0/a=5.191(16)`$) with the Symanzik improved gauge action and Asqtad quark action \[?\]. At $`\beta =6.0`$ we used for the thin link action data from ref. \[?\]. The quenched points show clearly the good scaling properties of the HYP action. The available dynamical data are in agreement with the quenched ones.

## 3 The New Algorithm

We propose a two–step updating to simulate the system described by an action $`S`$ as given in eq. (1)

* update a set of thin links $`\{U\}\{U^{}\}`$ in a way that satisfies detailed balance for $`S_g(U)`$
* accept/reject the proposed change with acceptance probability

$$P_{\mathrm{acc}}=\mathrm{min}\{1,\frac{det\left[Q^{}(V^{})Q(V^{})\right]}{det\left[Q^{}(V)Q(V)\right]}\},$$
(3)
where $`V^{}`$ are the smeared links constructed from the updated thin links.

Since the smeared links $`V`$ are constructed deterministically form the thin links $`U`$, they are not dynamical variables and an ergodic update of the $`U`$ links will correctly simulate the system. The proof that the above two–step updating satisfies detailed balance for the action $`S`$ follows closely the proof given in \[?\] for a similar two–step decomposition of the action.

In the practical implementation of the algorithm, our choice for $`S_g`$ is the Wilson plaquette action and the updating in step 1) is performed either with microcanonical overrelaxation \[?,?\] or with Cabibbo-Marinari heatbath \[?\]. We choose the subset of thin links to be updated differently for the overrelaxation and for the heatbath steps. For the overrelaxation we reflect all the links within some finite block of the lattice. The location of the block is chosen randomly but its dimensions are fixed. We choose with probability 1/2 a given sequence of reflections within this block or with equal probability the reversed sequence. The sequence has to be reversed with respect to the direction and location of the thin links and with respect to the index of the SU(2) subgroup used in the reflection step. For the heatbath we choose the subset of thin links randomly, i.e. we choose a random direction, a random parity, a set of random sites and a random sequence of SU(2) subgroups. The probability of generating a given sequence of thin link updates or the reversed sequence is then the same, both for overrelaxation and heatbath. From this it follows that detailed balance with respect to the pure gauge action $`S_g`$ is satisfied for both updatings. We chose to update a contiguous block of links with overrelaxation in order to propagate a change through (a portion of) the lattice. We observed that a random sequence of overrelaxed reflections is less efficient, especially for changing the topological charge of the configurations. This can be understood by thinking of instantons as extended objects. One has to change an entire region of the lattice to destroy or to create them.

In step 2) we use a stochastic estimator to evaluate the ratio of fermionic determinants

$`P_{\mathrm{acc}}^{}`$ $`=`$ $`\mathrm{min}\{1,\mathrm{exp}\mathrm{\Delta }S\}`$ (4)
$`\mathrm{\Delta }S`$ $`=`$ $`\xi ^{}[Q^{}(V^{})Q(V^{})Q^{}(V)Q(V)]\xi ,`$ (5)
where the vector $`\xi `$ is generated according to the probability distribution

$$P(\xi )\mathrm{exp}\{\xi ^{}Q^{}(V^{})Q(V^{})\xi \}.$$
(6)
The key feature is that smeared links constrain the statistical fluctuations of the stochastic estimator and make the algorithm effective, as we will show in the next section for the case of HYP staggered fermions.

Yet another interesting question is about the actual value of the ratio of fermionic determinants in eq. (3). This value is the natural limit on the performance of the algorithm, something like the Carnot efficiency of an heat engine. No matter how the accept/reject step 2) is performed, the detailed balance condition implies that the maximal acceptance is given by the Metropolis acceptance eq. (3) 111We thank U. Wolff for emphasizing this point.. We use a stochastic estimator to compute the acceptance and this reduces in principle the efficiency of the algorithm. The actual value of the ratio of fermionic determinants depends on the change in the gauge configuration. A larger number of changed thin links $`U`$ leads to a smaller value of the determinant ratio.

## 4 Four flavors of HYP staggered fermions

In the following we describe the application of the new algorithm to the case of $`n_f=4`$ flavors of staggered fermions coupled via hypercubic (HYP) smeared links. The matrix $`Q`$ in eq. (1) is given by

$`Q_{i,j}(V)`$ $`=`$ $`2(am)\delta _{i,j}+D_{i,j}(V),`$ (7)
$`D_{i,j}(V)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\eta _{i,\mu }(V_{i,\mu }\delta _{i,j\widehat{\mu }}V_{i\widehat{\mu },\mu }^{}\delta _{i,j+\widehat{\mu }}).`$
If we take the matrix $`Q^{}Q`$ to be defined on the even sites of the lattice only \[?\], eq. (1) describes four flavors of staggered fermions coupled via HYP smeared links $`V`$ ($`\eta _{i,\mu }`$ are the staggered phases). To further reduce the fluctuations of the stochastic estimator in the accept/reject step eq. (4), we remove from the fermion matrix its most ultraviolet part by decomposing it as

$`Q(V)`$ $`=`$ $`Q_r(V)A(V),`$ (8)
$`A(V)`$ $`=`$ $`\mathrm{exp}[\alpha _4D^4(V)+\alpha _2D^2(V)].`$ (9)
This way we achieve that an effective action

$$S_{\mathrm{eff}}=2\alpha _4\mathrm{Re}\mathrm{trD}^42\alpha _2\mathrm{Re}\mathrm{trD}^2$$
(10)
is removed from the fermionic determinant. The acceptance probability in eq. (4) becomes

$`\stackrel{~}{P}_{\mathrm{acc}}`$ $`=`$ $`\mathrm{min}\{1,\mathrm{exp}[S_{\mathrm{eff}}(V)S_{\mathrm{eff}}(V^{})]\mathrm{exp}\mathrm{\Delta }S_r\}`$ (11)
$`\mathrm{\Delta }S_r`$ $`=`$ $`\xi _r^{}[Q_r^{}(V^{})Q_r(V^{})Q_r^{}(V)Q_r(V)]\xi _r,`$ (12)
where the vector $`\xi _r`$ is generated according to the probability distribution

$$P(\xi _r)\mathrm{exp}\{\xi _r^{}Q_r^{}(V^{})Q_r(V^{})\xi _r\}.$$
(13)
In practice, we start by generating a random Gaussian vector $`R`$ from which the vectors

$`\mathrm{\Phi }^{}`$ $`=`$ $`Q^{}(V^{})R,`$ (14)
$`X^{}`$ $`=`$ $`[Q^{}(V^{})Q(V^{})]^1\mathrm{\Phi }^{}.`$ (15)
are constructed. The vector $`\xi _r`$ in eq. (13) is then given by $`\xi _r=A(V^{})X^{}`$. The two terms in eq. (12) are calculated with the formulae

$`\xi _r^{}Q_r^{}(V^{})Q_r(V^{})\xi _r`$ $`=`$ $`\mathrm{\Phi }^{}X^{},`$ (16)
$`\xi _r^{}Q_r^{}(V)Q_r(V)\xi _r`$ $`=`$ $`X^{}A(V^{})A(V)^1Q^{}(V)Q(V)A(V)^1A(V^{})X^{}.`$ (17)
All the vectors $`\mathrm{\Phi }^{}`$, $`X^{}`$, $`\xi _r`$ and the traces in eq. (10) are defined on the even sites of the lattice only.

The real parameters $`\alpha _4`$ and $`\alpha _2`$ can be optimized to maximize $`\stackrel{~}{P}_{\mathrm{acc}}`$, eq. (11). We keep the choice $`\alpha _4=0.006`$, $`\alpha _2=0.18`$ from our previous work \[?\] where the values were found to be optimal for an action with smeared links obtained by three levels of ordinary APE blocking \[?\]. The computation of the matrix $`A`$ in eq. (9) requires the evaluation of the exponential function which we do by truncating the exponential series after 15 terms.

In Monte Carlo simulations we separate the measurements of the observables by $`N_{\mathrm{OR}}`$ overrelaxation and $`N_{\mathrm{HB}}`$ heatbath two–step updatings, changing $`t_{\mathrm{OR}}`$ and $`t_{\mathrm{HB}}`$ links respectively. On $`8^3\times 24`$ lattices we simulated the HYP action with gauge coupling $`\beta =5.2`$ and bare quark mass $`am=0.1`$. An approximate matching of the lattice spacing ($`a0.17\mathrm{fm}`$) and of the Goldstone pion mass ($`m_\mathrm{G}r_0=2.0`$) can be achieved by simulating the thin link action with parameters $`\beta =5.2`$ and $`am=0.06`$. In an attempt to make a reliable comparison of the computational cost of the two simulations, we give in the last column of table 2 estimates for the number of $`Q^{}Q`$ multiplications required to get an independent configuration. These estimates are obtained by considering the autocorrelation time $`\tau _{\mathrm{int}}(p)`$ for the plaquette. By independent configurations we mean configurations obtained every $`2\tau _{\mathrm{int}}(p)`$ updating steps.

We found that measurements with the HYP action separated by $`N_{\mathrm{OR}}=160`$ overrelaxation and $`N_{\mathrm{HB}}=80`$ heatbath updatings have $`\tau _{\mathrm{int}}(p)=15(7)`$. Each overrelaxation updating changes $`t_{\mathrm{OR}}=128`$ and each heatbath updating changes $`t_{\mathrm{HB}}=200`$ links with an acceptance of $`20\%`$, corresponding to an effective change of 15% of all the links between measurements. Each updating, with overrelaxation or heatbath, requires one inversion of $`Q^{}(V)Q(V)`$, see eq. (15), performed with one conjugate gradient iteration (CG). Each CG iteration requires with these parameters on average 77 multiplications by $`Q^{}(V)Q(V)`$. In addition, each updating requires twice the computation of the matrix $`A(V)`$, see eq. (17), for a total of 60 multiplications by $`Q^{}(V)Q(V)`$. The recomputation of the HYP links once 200 dynamical thin links are changed has a cost on this lattice volume comparable to 33 multiplications by $`Q^{}(V)Q(V)`$. All together, we need 170 multiplications by $`Q^{}(V)Q(V)`$ for each updating.

We simulated the thin link action with HMC. We found that measurements separated by one HMC trajectory of unit time length have $`\tau _{\mathrm{int}}(p)=12(5)`$. One trajectory requires about 67 CG iterations, each of which needs with the given parameters on average 135 multiplications by $`Q^{}(U)Q(U)`$.

In summary, an independent configuration for the HYP simulation costs $`1.2(6)\times 10^6`$ multiplications by $`Q^{}(V)Q(V)`$ and for the thin link simulation $`0.22(8)\times 10^6`$ multiplications by $`Q^{}(U)Q(U)`$. This gives a factor 5(3) for the HYP computational cost over thin link HMC. To achieve the same improvement of flavor symmetry and scaling of the HYP action, the lattice spacing of thin link action simulations should be reduced at least by a factor 2. The overall gain of simulations with the HYP action is evident.

A drawback of our algorithm is that as the lattice volume increases the numbers of the updated links $`t_{\mathrm{OR}}`$ and $`t_{\mathrm{HB}}`$ have to be kept unchanged, consequently the numbers of updatings $`N_{\mathrm{OR}}`$ and $`N_{\mathrm{HB}}`$ have to be scaled with the volume to keep the autocorrelation times unchanged. Since each updating requires the evaluation of the fermionic action this gives a volume square dependence for the cost of the algorithm at fixed lattice spacing. On the other hand as the continuum limit is approached the numbers of links $`t_{\mathrm{OR}}`$ and $`t_{\mathrm{HB}}`$ which can be effectively updated scale: the physical volume of the updated region is constant. Considering our simulation of the HYP action at $`\beta =5.2`$ and $`am=0.1`$, we can update a physical volume $`0.027(1)\mathrm{fm}^4`$ with one overrelaxation updating and $`0.042(2)\mathrm{fm}^4`$ with one heatbath updating with 20% acceptance rate. Simulating the HYP action at $`\beta =5.0`$ and $`am=0.1`$, with a lattice spacing which is 30-40% larger, we can update a physical volume $`0.034(6)\mathrm{fm}^4`$ with one overrelaxation updating and $`0.056(10)\mathrm{fm}^4`$ with one heatbath updating with 20% acceptance rate.

To improve the efficiency of the HYP update one can improve the stochastic estimator needed in the acceptance step, adding a $`D^6`$ term in the definition of the matrix $`A`$ eq. (9). Furthermore, the pure gauge action used in the heatbath and overrelaxation updatings can be changed to make a better proposal for the acceptance step and the convergence of the exponential series calculation for the matrix $`A`$ can be made faster. Some of these improvements are being discussed in a forthcoming publication \[?\].

## 5 Conclusions

We presented a new simulation algorithm applicable to QCD actions where the fermions are coupled via smeared links. The algorithm is based on the standard overrelaxation and heatbath algorithms. We used it to simulate four flavors of staggered fermions with HYP smeared links. The cost of the algorithm has a volume square dependence at fixed lattice spacing due to the finite volume of the lattice region that can be updated, but the physical volume of the updated region remains constant as the continuum limit is approached. A comparison with the time cost of standard staggered action simulations with Hybrid Monte Carlo shows that the new algorithm is a factor 2–8 slower on lattices of physical volume of about $`10\mathrm{fm}^4`$ at $`m_Gr_0=2.0`$. The improvement in flavor symmetry and scaling of the HYP action allows to gain a factor 2 in the lattice spacing $`a`$. Given that the present cost of QCD simulations grows at least as $`1/a^6`$ there is a consistent improvement with the new algorithm.

An interesting point is that the algorithm could be used to simulate chiral fermionic actions like overlap or perfect actions. The presence of smeared links is the key feature which makes the algorithm efficient.

Acknowledgement. We are indebted to M. Hasenbusch for discussions that led to the new algorithm. We benefited from many discussions with Prof. T. DeGrand during the course of this work. We thank Prof. D. Toussaint for sharing with us the quenched Asqtad data of the MILC collaboration. The simulations were performed in scalar and parallel mode on the beowulf cluster of the high energy theory group and on the linux farm of the high energy experimental group at University of Colorado. We would like to thank the system administrator, D. Johnson, for his technical support. This work was supported by the U.S. Department of Energy. Finally we thank the MILC collaboration for the use of their computer code.