Improvement of N_f=3 lattice QCD with Wilson fermions and tree-level improved gauge action

We determine the parameter c_SW required for O(a)-improvement of the three flavor Wilson fermion action together with the tree-level Symanzik improved gauge action. The standard improvement condition is employed for a range of couplings. Additionally, we perform a check of the volume independence of c_SW and provide a preliminary estimate of the lattice spacing at our largest values of g_0^2.


Introduction
The continuum limit is an essential part of lattice QCD calculations. In order to control this extrapolation in the lattice spacing a, it is advisable to work with a discretization whose leading cut-off effects are O(a 2 ). The theory and reduction of scaling violations in on-shell quantities is well-established [1,2,3] and for Wilson fermions [4] O(a) improvement can be achieved by adding a single dimensionfive operator to the action [5]. This requires the non-perturbative tuning of its coefficient c SW and can be performed using a standard procedure [6,7]. Simulations with the resulting action in two-flavor flavor QCD have exhibited the expected moderate scaling violations [8].
Here we present a determination of the coefficient c SW using the standard procedure for N f = 3 flavor QCD with the tree-level improved Lüscher-Weisz gauge action [9]. We choose this gauge action because it has been demonstrated to possess superior scaling properties in pure gauge theory [10].
This paper is organized as follows. After defining our lattice regularization in Sec. 2, we summarize the standard improvement programme in Sec. 3. In Sec. 4 details of the numerical computations are given together with the resultant interpolating formula for the improvement coefficient c SW (g 2 0 ). We also give results for a renormalized quantity which suggests that the range of bare couplings covered by our simulations extends at least to a lattice spacing of a ≈ 0.09 fm.

Lattice setup
The O(a)-improved Wilson Dirac operator [4,5] is given by with ∇ µ and ∇ * µ the covariant forward and backward derivatives, respectively, and F µν the standard discretization of the field strength tensor [11]. The bare mass m 0 will be replaced below by the hopping parameter κ with m 0 = (κ −1 − 8)/2.
The gauge action S G contains sums over all oriented 1 × 1 plaquettes as well as all 1 × 2 rectangles, which are denoted by the sets S 0 and S 1 , respectively, where β = 6/g 2 0 , c 0 = 5/3 and c 1 = −1/12. The weight factor w k (C) is set to unity for all loops away from the boundaries. Schrödinger functional [12,13] boundary conditions are imposed on the gauge fields such that w 0 (C) = 1 2 , all links in C are on a boundary 1, otherwise (2.3) and all links in C are on a boundary 3 2 , C has exactly two links on a boundary 1, otherwise.
This corresponds to 'Choice B' of Ref. [14] and guarantees boundary O(a) improvement at tree-level of perturbation theory. We also set the fermionic boundary counter term according to the 1-loop formula [14] c Finally, the values of the spatial links at the boundaries are fixed to while the fermion fields satisfy periodic boundary conditions in the spatial directions and the standard Schrödinger functional boundary conditions in time.

Improvement condition
The standard O(a) improvement programme relies on the PCAC relation, which involves the improved axial-vector current (A I ) a µ and the pseudoscalar density P a given by In the unimproved theory the unrenormalized PCAC relation is violated by terms of O(a). By using three different choices of x, O and O ′ one can define c SW and c A , requiring that m is the same in all three cases. In this situation Eq. 3.1 holds up to O(a 2 ). This method has been applied to the N f = 0, 2, 3, 4 cases for a variety of different actions [6,7,15,16,17,18]. For technical reasons, we use lattices with T = 2L − a. However, this additional O(a) effect is irrelevant, as the determination of c SW using our improvement condition is ambiguous at O(a). We employ the PCAC relation with boundary operators O and O ′ on time slices x 0 = 0 and x 0 = T , respectively, The correlation functions which enter the PCAC relation Eq. 3.1 are then As has been suggested in Ref. [6], effective masses M(x 0 ) defined by correspond to a particular choice of c A in the improved currents. Since M(x 0 , y 0 ) renormalizes multiplicatively, it is a useful quantity to define the quark mass at fixed β and c SW . Specifically, we take as our definition of the quark mass Tuning M ≈ 0 defines the values of κ at which we impose the improvement condition.
As discussed above, this improvement condition requires the difference between two masses to vanish. In addition to M(x 0 , y 0 ), a second mass M ′ (x 0 , y 0 ) is considered where r and s in Eq. 3.7 are replaced by their primed counterparts. These masses are evaluated at x 0 = 3T /4 and y 0 = T /4. However, because we have T = 2L − a, we round the two arguments towards the center of the lattice.
The parameter c SW is chosen such that the difference between these two masses ∆M is equal to its tree-level value ∆M (0) The numerical value of ∆M (0) depends on the lattice geometry and can be computed using the solution to the Dirac equation as given in Sec. 6.2 of Ref. [19]. For the 15 × 8 3 lattices a∆M (0) = 0.000393. Additionally, this offset may be obtained from measurements on free gauge fields. A stringent test of our entire workflow is the reproduction of the ∆M (0) obtained from analytic calculations using simulations at large values of β.
In principle it would be preferable to keep the physical size of the system L constant as the continuum limit is approached. This, however, turns out to be very costly in practice, because autocorrelations associated with the topological charge sectors quickly become very large as the lattice spacing is lowered. We therefore opt to impose the improvement condition at fixed L/a, where the contribution of non-zero topological charge sectors decreases rapidly in the continuum limit. In Sec. 5.3 we will show that the effect of the finite volume does not seem to be relevant at the current level of accuracy.

Simulations
For the simulations we use the openQCD code, which is publicly availible online 1 and implements the lattice setup of Sec. 2 for several simulation algorithms. These algorithms are the subject of Ref. [20] so we restrict ourselves to a brief summary here.
We employ the HMC algorithm [21] with a twisted-mass Hasenbusch frequency splitting [22,23] for a doublet of two of the three degenerate quarks. For most ensembles with β ≤ 3.5 twisted mass reweighting [24] is used, i.e. we simulate with a small twisted mass µ = 0.001 and then include a stochastically estimated reweighting factor to correct for this in the measurement. This significantly increases the stability of the simulation. The third quark is simulated using the RHMC algorithm [25]. In all cases, we use a nine pole rational approximation in the interval [0.02, 7.2] and correct for the rational approximation with an additional stochastically estimated reweighting factor. All fermion determinants are factorized using even-odd preconditioning.
For the entire range of couplings, we generated 15 × 8 3 lattices at a variety of c SW values listed in Table 1. In addition, 23 × 12 3 ensembles at β = 3.8 have been generated as finite volume checks. To obtain a preliminary estimate of the lattice spacing at β = 3.3 and 3.4, we performed some L = T runs (using a modified version of the openQCD code) with L/a = 8 and zero boundary fields. These runs used both our discretization and the one of Ref. [26] and are summarized in Table 3.
For the L/a = 8 lattices at the four smallest values of β we use a three-level hierarchical integration scheme in which the outermost level employs the second order Omelyan-Mryglod-Folk (OMF) integrator, while the two inner levels use the fourth order OMF integrator [27]. The force from the pole closest to the origin in the rational approximation is integrated on the coarsest timescale, the remaining fermion forces on the intermediate timescale, and the gauge force on the finest timescale. Four or five outermost iterations together with a single iteration of the remaining two integrators typically achieve ≈ 90% acceptance for molecular dynamics trajectories of length τ = 2.
For the L/a = 12 runs at β = 3.8, the remaining L/a = 8 runs and the T = L runs, a two level scheme with both levels using fourth order OMF was found to be effective, with between 5 and 8 steps in the outermost integrator and a single inner iteration.
The integration schemes discussed above were very stable in all cases, resulting in small Hamiltonian violations. The most difficult simulations were those with L/a = 8, T = 2L − a at β = 3.3 at the smallest value of c SW , but even those had only about 0.6% trajectories with ∆H > 10. At c SW = 2.1 this falls already to 0.2% and we had 0.08% of such events at c SW = 2.4. At the smallest c SW = 1.7 for β = 3.4 these trajectories occurred with 0.1%, a percentage rapidly falling for larger values of c SW and β.

Results
We finally come to the results of our simulations. Our analysis strategy is discussed in Sec. 5.1 and results for c SW (g 2 0 ) are collected in Sec. 5.2, with a finite volume check in Sec. 5.3. A preliminary scale determination is given in Sec. 5

.4.
A summary of the ensembles generated for the determination of c SW appears in Tab. 1, where we also give the total statistics accumulated over several replica. These ensembles consist of L/a = 8 lattices used for our final result as well as those used for the finite volume check. The ensembles generated for the preliminary scale setting are discussed in Sec. 5  Values of the optimal c SW for which a∆M = a∆M (0) . Each value is from a linear interpolation of L/a = 8 data at a fixed value of β.

Analysis
For each value of the coupling, we have to find the value of c SW and the quark mass for which the improvement condition is satisfied. It has already been shown in the past that the condition M = 0 does not have to be achieved to a very high accuracy [7,17] and we therefore require |aM| < 0.015 for each individual point. We then measure a∆M for several choices of c SW and interpolate linearly to find the point with ∆M = ∆M (0) . The measurements of the required fermionic correlation functions (Eq. 3.5) are separated by one trajectory of length τ = 2. Using the methods and software of Ref. [28], we determine the integrated autocorrelation times of the observables aM and a∆M and find them to be at most around 8 units of molecular dynamics time such that all ensembles provide at least 500 independent measurements and thus a reliable determination of the errors. This is confirmed by the normal distribution of mean values on single replica.
Also, by studying the observables defined through the Wilson flow [29]the action density constructed from smoothed links and the topological chargewe ensure that field space is sufficiently sampled in all simulations. In particular the topological charge is known to cause problems in the continuum limit [30]. However, since we keep L/a fixed instead of L, the contribution of sectors of non-zero topological charge diminishes rapidly as a → 0.
The integrated autocorrelation times of the smoothed action and the topological charge at various smoothing ranges do not exceed 100 units of molecular dynamic time in any of our simulations. We therefore conclude that also for these observables configuration space is sampled sufficiently.

c SW for L/a = 8
The results for the improvement condition as a function of c SW as well as the resultant linear interpolations are shown in Fig. 1 for six values of β. The level of statistics is such that c SW is determined with better than 3% precision in all cases. The optimal c SW values at fixed β are collected in Tab. 2 and shown in  Table 2. We see that the fit describes the data well and that the data approaches the known 1-loop result at large β, which we use to constrain the interpolating curve.

Finite volume check
Although c SW has typically been determined at fixed L/a = 8, it is interesting to assess the impact of this choice on the final result. As stated above, it would actually be preferable to keep L fixed when imposing the improvement condition. Smaller volumes are advantageous as they are computationally cheaper and have a larger slope in ∆M vs. c SW , resulting in a clearer signal. Also, the problem of autocorrelations associated with freezing topological modes is much reduced. However, small volumes give rise to potentially large O(a) effects in c SW . We detail here a single check at β = 3.8 between L/a = 8 and L/a = 12. The result is shown in Fig. 3. Within the statistical accuracies, the two volumes give the same value of c SW . With c SW = 1.64(2) for L/a = 8 and 1.61(5) for L/a = 12, To quantify finite volume effects, a value from simulations at 23 × 12 3 is given, with the value for g 2 0 slightly shifted to the left for clarity.
one can conclude that at least up to that lattice spacing, the smaller lattices are a good choice for the determination of c SW . The point at L/a = 12 also does not deviate by more than one standard deviation from the interpolating curve plotted in Fig. 2. Therefore, given our statistical accuracy, the systematic effects from fixing L/a = 8 do not appear to be significant.

Preliminary scale determination
To estimate the lattice spacing we compute a quantity with our discretization as well as the discretization used by PACS-CS, where the scale is known from large volume simulations using the mass of the Ω baryon [31]. For comparison we use the coupling defined in Ref. [32], except with periodic spatial boundary conditions for the fermion fields, i.e. with θ = 0. This coupling is a renormalized quantity and at M = 0 depends only on L, up to scaling violations. This lattice size L at a given value of the coupling serves as the dimensionful quantity to set the physical scale.
Our results for this coupling are tabulated in Table 3. In both discretizations we expect cutoff effects of order O(ag 2 0 ) as boundary improvement for the gauge fields is implemented at tree level only. Furthermore, boundary improvement for the fermion fields is implemented only at 1-loop, resulting in additional O(ag 4 0 ) effects. We see that the β = 3.3 result for the coupling in our discretization lies above the value at a = 0.09fm, while the β = 3.4 is below. This suggests that the β corresponding to a = 0.09fm in our discretization is in the range [3.3, 3.4].

Conclusions
In this paper we have determined c SW for N f = 3 lattice QCD with the treelevel Symanzik-improved gauge action. The result of our determination is the interpolation formula which may be taken as a definition of the lattice action. While our determination was performed at fixed L/a = 8, we performed a finite volume check at β = 3.8 and found no significant change in c SW . In addition to a determination of c SW , we have calculated a renormalized Ldependent coupling at β = 3.3 and 3.4 and compared with an alternative N f = 3 discretization. This determination suggests that the bare coupling corresponding to a ≈ 0.09fm in our discretization is located in the interval β ∈ [3.3, 3.4], indicating that our c SW determination spans the range of desired lattice spacings.  Table 3: N f = 3 results for the Wilson flow couplingḡ 2 (L) from Ref. [32] using the Iwasaki gauge action and our discretization. We take the parameters and the scale determination from Ref. [31], while κ was taken from Ref. [20] and found to give aM ≈ 0. For these simulations only, we set T = L with boundary fields C k = C ′ k = 0.