Taming lattice artifacts with Pauli--Villars fields

As fermions are added to a lattice gauge theory, one is driven to stronger bare coupling in order to maintain the same renormalized coupling. Stronger bare couplings are usually associated with larger gauge fluctuations, leading to larger cutoff effects and more expensive simulations. In theories with many light fermions, sometimes the desired physical region cannot be reached before encountering a phase boundary. We show that these undesired effects can be reduced by adding Pauli--Villars fields. We reach significantly larger renormalized couplings while at the same time damping short-distance fluctuations of the gauge field. This may allow for controlled continuum extrapolations from large lattice spacings.


INTRODUCTION
There is wide interest in the non-perturbative properties of asymptotically free gauge theories at strong renormalized coupling. Lattice gauge theory simulations aim to investigate strong-coupling physics at large distances while starting at weak coupling at short distances. Tuning the bare coupling to zero decreases the lattice spacing and reduces cutoff effects but requires increasing the lattice volume to keep the physical volume sufficiently large. Various improvement programs, both perturbative and non-perturbative, attempt to reduce lattice artifacts so that simulations even at reasonable lattice volumes can predict continuum physics with high precision. This approach works well when the gauge coupling runs fast, as in QCD with a few flavors of light fermions.
As is well known, fermion fields screen the gauge interaction and thus slow down the running of the gauge coupling. Sufficiently many fermions even lead to the emergence of an infrared fixed point. In slowly running systems, computationally acceptable lattice volumes do not allow significant change from the bare to the renormalized coupling. In order to reach strong renormalized couplings, simulations must start with a large bare coupling. This, in turn, brings in large cutoff effects. Often the large ultraviolet fluctuations induce unwanted lattice artifacts like bulk phase transitions.
We illustrate the problem with Fig. 1. The top panel shows g 2 GF , the gradient-flow coupling 1 at fixed lattice distance as a function of the bare coupling β for SU (3) gauge theory coupled to nHYP staggered fermions with N f = 0, 4, 8, and 12 (continuum) flavors. Each theory except the pure Yang-Mills theory exhibits a first-order bulk transition to an S4 phase 2 at strong coupling. The leftmost point in each dataset is slightly to the right of the phase transition for that model. One can see that the largest g 2 GF that can be reached before encountering the phase transition gets smaller with increasing N f . The bottom panel of Fig. 1 shows the corresponding plaquette values (normalized such that the maximal value is 3). Comparing the two panels reveals that, as the number of fermions increases, reaching the same value of g 2 GF (if possible) is accompanied by a smaller plaquette value. In other words, it comes at the price of stronger shortdistance gauge field fluctuations.
Small plaquette values typically signal a difficult simulation. Worse is the presence of the phase boundary. There is the possibility that simulations intended to study the properties of the Gaussian fixed point at g 2 0 = 0 but performed near the first-order transition will reflect the properties of the bulk transition. In any case, bulk phase transitions limit the physical range that lattice simulations can probe.
Cutoff effects will be reduced if we can suppress ultraviolet fluctuations without affecting the infrared physics. In this paper we show that adding to the lattice action a set of heavy Pauli-Villars (PV) fields achieves this goal. These are bosonic Dirac fields that we choose to have the same lattice action as the fermions, but with heavy masses fixed at am PV = O(1). These PV fields play a role similar to that of the continuum PV fields in regularizing fermions at large momenta [5]. While fermions screen the gauge coupling, the PV fields anti-screen. At short distances r < ∼ 1/m PV they will enhance the running of the coupling. At large distances r a the running of the gauge coupling is physical, and will not be altered by the presence of the PV fields.
Matter fields whose mass is fixed in lattice units, whether fermions or bosons, decouple in the continuum limit where the lattice spacing a → 0 or the lattice correlation length ξ → ∞. The effect of the PV bosons can be summarized as an addition to the lattice action. They generate a local effective gauge action [6] that can balance the screening effect of the fermions at short distances. This recalls the old idea of gauge fields wholly generated by heavy matter fields [7,8], but here the goal is only to compensate for the short-distance effects of lattice fermions. While continuum PV fields are introduced in one-to-one correspondence with the fermions in order to cancel UV divergences, on the lattice we may allow the number of PV fields to be a free parameter. The induced term in the gauge action is simply proportional to the number of PV fields. 3 This paper is organized as follows. In Sec. II we present the PV fields that accompany the staggered fermions in our simulations. We also review the idea that massive fields can be considered as generating an effective action for the gauge field. Section III gives an account of our experimentation with the PV fields, and in Sec. IV we give our conclusions.

A. Lattice definitions
In this paper we use a gauge action that includes fundamental and adjoint plaquette terms with couplings β ≡ β F and β A , respectively, related by β A /β F = −0.25 [1]. The gauge links in the staggered operator D are nHYP-smeared links [13,14] with smearing parameters α = (0.5, 0.5, 0.4). This lattice action has been used in several studies of 8-and 12-flavor systems, including explorations of the phase diagram [1,2], the discrete β function [15,16], the scale-dependent mass anomalous dimension [17], and the large-scale spectrum studies of the LSD collaboration [18,19].
The finite volume gradient-flow coupling is defined by [20][21][22][23] where N c = 3 and the term 1/(1 + δ) corrects for the gauge zero modes due to periodic boundary conditions [23]. We evaluate g 2 GF at flow time fixed by √ 8t = cL; throughout this paper we use c = 0.45.
In Eq. (2.1), E(t) is the (flowed) energy density. Different operators and flows, like the plaquette and clover operators, or Wilson and Symanzik flows, give slightly different g 2 GF , but for √ 8t = 0.45 L this variation is usually mild. In the following we define g 2 GF using Wilson flow and the clover operator. Since we use the same flow and operator throughout this work, the specific choice should not affect the qualitative results.
We define the staggered lattice Dirac operator as where D is the single-component, massless, antihermitian staggered operator. The fermion action constructed using this operator gives rise to N f = 4 continuum flavors when the bare gauge coupling is tuned to the Gaussian fixed point g 2 0 = 0 (equivalently, β → ∞), provided that the fermion mass m is accordingly tuned to zero in lattice units.
The squared determinant of M can be represented using a single pseudofermion field φ, 3) The path integral that includes (DetM ) 2 gives N f = 8 flavors in the continuum limit. By restricting M † M to the even sites only, one obtains N f = 4 flavors in the continuum limit.
Pauli-Villars fields have the opposite statistics. The bosonic inverse determinant is represented similarly as where we define Unlike the pseudofermion action, the PV action is local. As above, we restrict M † PV M PV to the even sites only. Including N PV copies of such PV fields in Eq. 2.4 would amount to 4N PV bosonic Dirac fields in the continuum limit if we were to tune the PV mass to zero. We stress, however, that the PV mass will always be held fixed in lattice units, in order to decouple the PV fields in the continuum limit.

B. Integrating out heavy matter fields
Matter fields with mass of the order of the cutoff decouple from the infrared dynamics, but they generate an effective gauge action, affecting the bare gauge coupling. N s staggered fermions give an induced contribution to the gauge action (2.6) The full effective action for the gauge field is thus For thin (unsmeared) links, an expansion in powers of 1/m gives where the sum over C is a sum over all closed loops of length originating from lattice site x, and E C = 1 4 Tr (γ µ1 · · · γ µ ) is a sign factor that depends on the geometry of the loop. For heavy staggered fermions interacting with thin (unsmeared) gauge links the leading term in S ind is a plaquette term with coefficient (2.9) The coefficient of F 2 µν from all the loops (representing the leading term in an expansion in the lattice spacing) can be expressed as a new inverse bare gauge coupling β eff = β + β ind . Smeared link actions generate a smeared plaquette at leading order that also corresponds to an induced action characterized by a shift in β.
Ref. [24] investigated the finite temperature phase transition of various systems (with thin links) and concluded that fermions with am > ∼ 0.1 are well described by an effective gauge action. The possibility of replacing the entire gauge action with an effective action generated by fermions was explored in [7], while [8] considered heavy scalar fields in the adjoint representation.
If we replace the fermions by PV bosons, the only change is that the overall sign on the right-hand side of Eq. (2.8) is flipped. For equal masses m PV = m, the PV bosons exactly cancel out the S ind generated by N s = N PV staggered fields. Similarly, in continuum perturbation theory, and for scales µ m PV = m, the contribution of the PV bosons to the one-and two-loop beta function would be exactly opposite to that of the fermions.
In the next section, we turn to a numerical study of the effect of PV bosons with masses am PV = O(1).

III. NUMERICAL SIMULATIONS WITH PAULI-VILLARS FIELDS
We have explored the effect of the PV fields in the SU(3) gauge theory with N f = 12 massless fermions, using a staggered action with nHYP smearing (see Sec. II A). This theory exhibits a first-order transition into the S4 phase [1][2][3][4], which persists for all the combinations of PV fields that we have tried out.
In Fig. 2 we show the effect of a growing number of PV fields, always keeping am PV = 0.5. We consider N PV = 0, 2, 4, 6, and 8. The lattice volume is 8 4 . The top panel of Fig. 2 shows g 2 GF as a function of the inverse bare gauge coupling β. A striking feature is how much fixed physics, defined by a given value of g 2 GF , gets pushed back towards larger β. This is thanks to the anti-screening effect of the PV fields. For example, let us consider g 2 GF = 11. Without PV fields, this value is achieved for β 2.8 (see Table 1). With N PV = 2, am PV = 0.5, the same renormalized coupling is attained for β 6.8. With each further addition of 2 to N PV , the value of β shifts by about 2.2. This behavior is qualitatively consistent with the predictions of the hopping-parameter expansion, see Sec. II B.
The bottom panel of Fig. 2 shows that increasing the number of PV fields while keeping the renormalized coupling fixed is accompanied by increasing plaquette values. This demonstrates that the short-distance fluctuations of the gauge field are damped. The system without PV bosons has plaquette expectation value around 0.6 (out of 3) just above the bulk transition. Such rough gauge fields would not be considered acceptable in QCD simulations.
Another striking feature seen in the top panel of Fig. 2 is the larger range of renormalized gauge couplings made accessible when PV fields are added. As in Fig. 1, the  (3) theory with N f = 12 fermions in the chiral limit. All systems exhibit a first-order transition to the S 4 phase. As before, the leftmost point in each set is a little to the right of the phase transition; data points inside the S4 phase are not shown.  leftmost point in each data set is just to the right of the phase transition into the S4 phase. As we increase N PV from left to right, the corresponding maximal values of g 2 GF are approximately 14, 17, 20, 25, and 29. 4 The largest renormalized coupling attainable with N PV = 8 is about twice that without any PV fields. Fig. 3 is similar to Fig. 2, but this time we keep the number of PV fields fixed at N PV = 4 and we consider three values of the PV mass, from large to small: am PV = 1.0, 0.75, and 0.5. As expected, the effect of the PV fields increases with decreasing mass. Comparing Figs. 3 and 2 we deduce, for example, that N PV = 4 with mass am PV = 0.75 has roughly the effect of N PV = 3 with m PV = 0.5. If parameters of the simulation require it, one can always raise the PV mass to am PV = 1, and then increase N PV till the desired anti-screening effect is achieved. (There is no reason to consider heavier PV bosons, since am PV = 1 already sets the PV mass equal to the cutoff.) We stress again that the PV mass should be held fixed in lattice units when taking the continuum limit.

NPV
The reader may wonder if PV bosons with mass am PV = 0.5 are heavy enough to decouple. First, we have found that the mass of the PV pion-the ghost pion made of two PV bosons-is roughly equal to 2m PV . This implies that the couplings between remote sites in S ind decrease exponentially with a decay rate ∼ 2m PV .
We can also examine this question from the point of view of the infrared scales. In systems where chiral symmetry is not broken, the only dimensionful parameters, besides m P V , are the fermion masses (denoted generically as m) and the lattice volume. Fermion bound states are expected to scale with m, so as long as m P V m and m P V 1/L, the PV bosons decouple. This applies both in mass-deformed conformal systems and in QCD-like systems at high temperature, where there is no condensate. In particular, if simulations are performed in the chiral limit, the only condition is m P V L 1. Where chiral symmetry is broken, an additional condition m P V Λ QCD has to be imposed. In practice one can check that m P V M H where M H is the heaviest physically accessible hadron of the system. The flow-time evolution of g 2 GF (t; L) is always characterized by three regions. It starts with an initial rise at t/a 2 < ∼ O(1) that depends strongly on the specific action, flow, and operator choice. This is followed by a physical region, and finally the flow saturates due to the finite volume. In Fig. 4 we show the t dependence of g 2 GF (t; L) for β just above the S4 transition, as shown in Fig. 3. In the absence of any PV fields the "knee" where g 2 GF turns from the initial short-distance rise to the physical (and, here, decreasing) part occurs at t/a 2 0.6. The inclusion of PV fields induces a short-distance gauge interaction that affects the small-t behavior. In the flows with 4 PV fields the turning point between the two regions moves to the right as m PV is lowered, while the maximum of g 2 GF grows. For c = 0.45 and L = 8 we read off the gradient-flow coupling at t/a 2 1.62. This small volume obviously does not leave much room for a physical region. In Fig. 4 we also show the flows for L = 12. Here the gradientflow coupling is determined at t/a 2 3.65, which is now well beyond the turning point of the flows for all values of m PV . Both the maximum of g 2 GF (t; L) and its value at √ 8t = 0.45 L do not change much with the volume. 5 Finally, we turn to the cost of the PV fields. As long as the number of PV bosons is not extremely large, adding them to the simulation is inexpensive. The PV action (2.4) requires only one inversion (with large mass) at the beginning of an HMC trajectory to generate each PV field Φ. During the molecular dynamics evolution of the gauge field, calculating the force due to the PV action requires mutiplication with the Dirac operator, but no inversion.
As a matter of fact, the PV fields can even make the simulations less expensive. In Table 1 we report the performance of simulations with different numbers of PV fields with m PV = 0.5. 6 As N PV is increased from zero, we change β to keep g 2 GF 11 fixed (cf. Fig. 2). We observe decreases in (1) the number of conjugate-gradient (CG) iterations and in (2) the average error |δH| in the molecular dynamics (MD) energy due to large forces in the equations of motion. Sharp decreases are seen in going from N PV = 0 to 2 despite an increased step size dτ = 1/N step . 7 A decrease in |δH| brings about improved acceptance and should allow yet larger values of dτ , further decreasing the computational cost.

IV. DISCUSSION
In this paper we investigated the possibility of countering the effective action generated by fermions at short distances by including heavy Pauli-Villars fields. The effective action generated by PV bosons has the opposite sign to that generated by fermions; while fermions screen the gauge coupling, PV bosons anti-screen. As long as their mass am PV is kept at O(1), the PV bosons decouple from the infrared dynamics.
Pauli-Villars fields are part of standard domain-wall fermion actions. Indeed, domain-wall fermions are generally known for their small cutoff effects. In this paper, we considered the use of a more general set of PV fields with staggered fermions. Similarly, one can add heavy PV bosons to a theory with Wilson fermions. We would always use the same lattice Dirac operator for the unphysical PV bosons as for the fermions (except of course with a heavy mass m PV ), anticipating that this is the 5 If, as in Fig. 2, we hold m PV fixed and vary N PV we observe a similar trend: again the turning point moves to the right with increasing N PV and the maximum of g 2 GF grows. 6 Chiral symmetry breaking was not observed in any of our data sets. This allowed us to work directly at am = 0. 7 Adding more than 4 PV fields seems to have no further effect. optimal way to reduce the cutoff effects of the fermions. One can also try using a different PV action; the PV fields would simply induce a different gauge action.
Increasing the number of heavy PV bosons (even beyond the number of fermions) enlarges the induced gauge action. As a consequence, the same physical regime is reached at smaller bare gauge coupling (i.e., larger β) where the gauge fields are smoother. The ability to work at small bare coupling while maintaining the desired renormalized coupling in the infrared means that cutoff effects are reduced. In a way, introducing many PV bosons brings the system closer to a perfect action.
Heavy PV bosons have an additional positive effect on systems that exhibit bulk first order phase transitions in strong coupling. These phase transitions are typically caused by large ultraviolet fluctuations and are not physical. Nevertheless they limit the range of possible renormalized couplings. While the additional PV bosons in our test study with N f = 12 flavors did not remove the bulk transition, they shifted the simulations to smoother gauge fields and opened up the range of accessible couplings. In our experiments we found about a factor of two increase in the maximal g 2 GF . Adding several PV fields to the simulation is inexpensive. In fact, the PV fields made our simulations less expensive: Since the gauge fields are smoother, the CG inversions needed for the fermion force converged significantly faster. Moreover, the smaller MD forces resulted in improved acceptance even with a larger step size.
The introduction of PV bosons can be particularly useful for slowly running systems where, for the renormalized coupling at the infrared scale to be large, normally the bare coupling would have to be taken large, too. PV bosons allow for independent control over the magnitude of discretization effects in such systems.
Our aim in this paper has been to show that the addi-tion of PV fields smooths the gauge field and allows the gradient flow to reach very strong couplings. We have demonstrated these claims using N f = 12 flavor staggered fermions. Any conclusions to be drawn from the data presented will, as usual, depend on the results of taking large-volume and continuum limits. The question of whether the N f = 12 theory is ultimately confining [25][26][27][28] or conformal [16,29] has been the subject of extensive study. If the downward trend of the flows in the physical region (Fig. 4) survives all relevant limits, this would signal the presence of an infrared fixed point at some smaller value of g 2 GF . PV bosons can find their use in QCD, too. For given computer resources, the only way to increase the physical volume is by increasing the lattice spacing. This comes at the price of growing discretization errors. Including PV bosons might allow for simulations with larger lattice spacings and hence with larger physical volumes, while still keeping the discretization effects under control. At the same time simulations with PV bosons could be computationally faster as well.