Setting the scale for nHYP fermions with the L\"uscher-Weisz gauge action

Lattice QCD calculations using gauge smearing for fermion kernels are computationally efficient. Hypercubic blocking (nHYP smearing) has been shown to reduce scaling errors. In this work we use an improved action for $N_f=2$ QCD, based on the L\"uscher-Weisz gauge action and clover-improved Wilson fermions with nHYP smeared gauge links. We perform a parameter scan in the region with lattice spacing between $0.066 \mathop{\hbox{fm}}$ and $0.115 \mathop{\hbox{fm}}$ and pion mass between $207 \mathop{\hbox{MeV}}$ and $834 \mathop{\hbox{MeV}}$. We determine the lattice spacing and pion mass as a function of the bare coupling parameters ($\beta$ and $\kappa$). The results are obtained from twenty-two ensembles on a $24^3\times 48$ lattice to percent level in statistical accuracy. The finite-volume effects for these ensemble are at the sub-percent level. From these measurements we produce easy-to-use parameterizations to help tune simulations with this action. The lattice spacing is fixed using a mass-independent procedure, by matching observables in the chiral limit. We also provide a parameterization for the chiral extrapolation which is universal and should hold for all discretizations of $N_f=2$ QCD.


I. INTRODUCTION
Quantum Chromodynamics (QCD) is the fundamental theory of the strong interaction, with only a few parameters: the strong coupling constant (g) for gluon-gluon and gluon-quark interactions, and the quark masses (m q ) (one for each quark flavor). In the energy region relevant for hadronic interactions the effective coupling is strong. Lattice QCD can be used to calculate QCD properties in this region via numerical simulations on a spacetime lattice with spacing a. Lattice QCD allows us to extract physical quantities in lattice units (or in dimensionless ratios) from Euclidean correlation functions. To relate such dimensionless quantities to physical quantities in physical units we must determine the lattice spacing. In practice, this is done by fixing one physical observable to its experimental value. Once the lattice spacing is fixed, the values of other observables can be converted to physical predictions. For this reason, scale setting is an important step in lattice QCD calculations.
Lattice QCD calculation are always done at finite lattice spacing. All predictions need to be extrapolated to the continuum limit, to eliminate the discretization errors. Generically discretizations errors are expected to vanish linearly with the lattice spacing as we take the continuum limit, but improved actions can be used to accelerate the convergence rate. The largest discretization errors are usually generated by the quark contribution to the action. Smearing the gauge fields was found to improve the scaling of hadron spectrum significantly [1][2][3][4][5]. With the introduction of analytical smearing [6], efficient dynamical simulations were made possible using Hybrid Monte Carlo (HMC) method [7]. Hypercubic blocking [2] was In this work, we consider an action for QCD with two mass-degenerate quark flavors using nHYP clover fermion discretization [8] for quark fields and Lüscher-Weisz gauge action [20]. We study the scale setting for this action to facilitate its wider usage. The plan of the paper is the following: In Sec. II, we give a largely self-contained description of the action. In Sec. III we present a brief review of the methodology for scale setting, our numerical simulation results, and finite volume effects. Smooth parametrization functions are extracted in Sec. IV, before conclusions in Sec. V.

II. ACTION
The lattice QCD action we use in this study is a function of the gauge links U µ and quark fields ψ, Note that the gauge part of the action S LW (U ) is a function of the thin links, whereas the links V that enter the fermionic matrix D clover are nHYP smeared links V (U ). In this section we review the details for each of the term in this action. operators beyond the usual plaquette, where we use the shorthand notation (3) The summation indices run in both positive and negative directions {±1, ±2, ±3, ±4}. For negative indices we have U µ (x) ≡ U −µ (x +μ) † . The set of indices in the parallelogram term has 16 elements, The three terms correspond to three types of gaugeinvariant lattice loops: square, rectangle, and parallelogram. These terms can be depicted graphically as, The coefficients in the action are taken from one-loop, tadpole-improved perturbation theory [21] where the measured expectation value of the plaquette is used to determine both the value of the mean link u 0 and the QCD coupling constant α s (the value of α s is around 0.3 in our simulations). This action has excellent rotational invariance properties in the static quark potential even on coarse lattices [21]. The scale for the pure gauge Lüscher-Weisz action has been set in an earlier study [22]. The fermion action is based on the clover discretization. The Dirac operator can be written as where D W is the Wilson term We work with r = 1. For the clover term, σ µν = 1 2 [γ µ , γ ν ] and the clover field F is constructed from the antihermitian and traceless part of the average of four plaquetes, where is represented pictorially in Fig. 1. The clover coefficient is set to its tree-level value c SW = 1.0 since we are using nHYP smeared links in the fermionic kernel, and non-perturbative determination of the clover coefficient found the corrections to be insignificant [9].
We now turn to the smearing function: nHYP smearing [8] refers to normalized hypercubic smearing which involves three consecutive levels of APE-style smearing [1], Here the conventional smearing notation using U n,µ to represent U µ (x) is adopted. The intermediate fieldsṼ andV are constructed such that the contributions to V n,µ are restricted to the thin links U appearing as edges in the hybercubes attached to the link V n,µ . The indices after the semicolons indicate the directions excluded from the sums. The operator P is a projection to U (3), The projection is non-singular by construction since (A † A) −1/2 is Hermitian and positive definite. The three U (3) projections render the nHYP smeared configurations very smooth, while keeping the smearing within a hypercube ensures that even short distance properties of the configurations are minimally distorted. Consequently, the fermion action is ultra local. Unlike HYP links which project to SU (3) rather than U (3), the nHYP links are differentiable, which makes them amenable to HMC algorithms. We use the standard HYP values α 1 = 0.75, α 2 = 0.6, and α 3 = 0.3 that have been tuned to minimize the plaquette fluctuations [2]. The idea of using smeared links in the fermion action to improve its scaling properties has been validated in a number of studies. In Ref. [3] it is shown that staggered fermions with HYP smearing has very good scaling properties, even better than other improved actions. In Ref. [8] and Ref. [9], the same nHYP fermion plus Lüscher-Weisz gauge action is explored on small lattices. The clover coefficient c SW is computed non-perturbatively and found to be close to the tree-level value of 1.0 with small corrections. In Ref. [4], tree-level clover fermion with stout links and tree-level Symanzink gauge action are used. A wide scaling region up to lattice spacing of 0.16 fm is found. Stout smearing [6] is an alternative projection method that is also differentiable. It can be used to replace the projection P discussed above, leading to the so-called HEX smearing. A two-level HEX smeared action is used in Refs. [4,5] to study continuum limit for hadron masses and quark mass renormalization.

A. Simulation parameters
We carried out dynamical simulations on 24 3 × 48 lattices with 6 values of β and up to 5 values of κ for each β, totaling 22 different ensembles. The gauge fields are updated using the HMC algorithm [7]. We employ the Hasenbush-Jansen multi-mass method [23] to speed up the calculation: by using two masses and adjusting the integration time steps for gauge and fermion force updates to keep the contribution balanced between the terms. We set the trajectory length to 1 and varied the number of integration steps to keep the acceptance rate above 90%. For each ensemble we generated 500 configurations using 2500 trajectories and saving every fifth trajectory. We then dropped the first 100 configurations for thermalization and kept the last 400 for calculations.
The list of all the ensembles is given in Table I along with other quantities determined in this work. The values for β were chosen to scan finely the 0.05 fm a 0.1 fm range. Since our simulations were run on a fixed sized lattice, to keep the finite volume errors under control, we used κ values such that 0.15 am π 0.3. For each ensemble we compute the pion mass using two point correlation functions with several sources per configuration, evenly spread out in the Euclidean time direction. The quark propagators were computed using our optimized GPU inverters [24]. We perform single exponential fits to extract the pion mass. The fit range for the pion correlators were determined by looking for plateaus in the effective mass plots for the correlators.
Since the observables are determined to percent, or better, precision level, we ensure that the systematic errors remain below percent level. In particular, we must assess finite volume effects for the pion mass. The finite volume effects for t 0 and w 0 , the other observables used in this study, are smaller [25].
Finite volume corrections to the pion mass can be related to the pion forward scattering amplitude [26]. Alternatively, such corrections can be computed using chiral perturbation theory in finite volume [27][28][29]. The two approaches were compared in Refs. [30,31], where the subleading effect in both approaches, the next leading exponential correction for Lüscher's approach and NLO effects for the χPT approach, were found to produce large corrections. However, as the χPT NNLO effects were found to be small, finite-volume χPT is expected to be reliable. To improve the convergence of Lüscher's method [26], a resummation method was proposed in Ref. [32] for finite-volume correction on the pion mass, where m(n) are multiplicities whose values can be found in Table 1 of Ref. [32]. The resummed correction shows very good agreement with the finite volume χPT expansion: when using the LO(NLO) χPT expression in the forward amplitude F π (iỹ), the results agree with 1-loop(2-loop) finite-volume χPT expansion. The resummed method was checked by direct 2-loop calculation of the finite-volume χPT and found to be accurate for m π L 2 [33,34]. Our values for m π L from Table I lie between 2.6 to 6.7. We use the χPT results to estimate the expected finite volume corrections: we interpolated the results from Table  3 of Ref. [32] to get the magnitude of these corrections for our ensembles. The resulting estimates are summarized in Table I. For our ensembles these corrections are at the level of 1% or below.

B. Lattice spacing methodology
There are several well-established methods used to determine the lattice spacing: the Sommer parameter [35], hadron masses [36], etc. While scale setting using hadron masses is conceptually straightforward, such quantities are often not determined with very high precision, and are computationally quite expensive. Additionally, systematic uncertainties associated with extracting hadron masses, such as excited state contamination and finite size effects, are not always easy to estimate. The Sommer scale r 0 is based on the calculation of the static potential from the gauge fields, which does not require the computation of expensive quark propagators. However, care must be taken to select the fit form and ranges. In addition to using observables which may be determined with good accuracy at minimal computational cost, it is often advantageous to use scales with a mild quark mass dependence [37]. The Wilson flow method proposed by Lüscher [38] satisfies all these requirements. This uses new parameters t 0 and w 0 [25,39], defined below, to set the scale.
The Wilson flow is essentially a smearing of the original gauge fields U µ (x) controlled by a flow time parameter t, not to be confused with x 4 , the Euclidean time coordinate. This method has the advantage that it is very straightforward to implement. The Wilson flow for gauge fields U µ (x) is defined by the first-order differential equation where S W is the standard Wilson gauge action (the first term in Eq. 2) evaluated using the smeared links V µ (x, t).
The smearing radius is proportional to √ t. We numerically integrate this equation for a flow time t 0 such that the smearing radius achieves a particular physical value.
To determine this point we monitor the energy density given by For the gluon field strength tensor F µν we use the clover discretization in Eq. 9 constructed from the smeared fields V µ (x, t). The parameters t 0 and w 0 are defined implicitly by the dimensionless, renormalized quantities An advantage of using t 0 and w 0 for scale setting is that, despite their relatively large auto-correlations, their stochastic errors are small and can be computed precisely and are computationally inexpensive. As t 0 and w 0 are not experimentally measurable quantities, in order to assign physical units we must compare to other lattice calculations in which a physical quantity has been used to set the scale (e.g. using hadron mass ratios). Hence, these Wilson flow observables offer an efficient way of determining the relative scale/lattice spacing between many simulations with differing gauge couplings. For N f = 2 we refer to Ref. [39] where t 0 and w 0 are reported in physical units at the chiral point. For a recent application of Wilson flow to set the scale for two-color QCD, see Ref. [40].

C. Lattice spacing results
We measure w 0 /a from Eq. 16 on multiple ensembles with the same gauge coupling β and extrapolate to the chiral point (m π = 0). We fix the scale by the value of w 0 at the chiral point w chiral 0 . We then use these measurements to find a parameterization of the lattice spacing w chiral 0 /a as a function of β. We are able to determine w 0 /a with very high accuracy at a few parts per-mille level.
To convert the lattice spacing results to physical units we use the N f = 2 value of w chiral 0 = 0.1776(13) fm from Ref. [39]. This value is extracted by computing the kaon decay constant in N f = 2 simulations [41]. Note that in this work, when we report the lattice spacing results in physical units, the bulk of the uncertainty in the lattice spacing comes from the uncertainty in w chiral 0 . The extrapolation of w 0 /a to the chiral point requires a detailed discussion. Following Ref. [39], we compute w 0 /a as a function of the variable y = (m π w 0 ) 2 using the ensemble data 1 . If the values of a are properly determined we can compute then w 0 as a function of y and these data points are expected to lie on a universal curve, up to discretization errors. In Ref. [39] this curve was found to be well described by a line, but our data-points have a wider range of y values and we see a deviation from the linear behavior. To account for the non-linear behavior, we perform a quadratic fit to w 0 /w chiral 0 as a function of y for each β and use the resulting function to interpolate w 0 /a as a function of y. 1 Notice that our definition of y here is different from the definition in Ref. [39] which is y = m 2 π t 0 .
To obtain the universal curve for w 0 (y) we use the following procedure. For each β we introduce a fit parameter (a/w chiral 0 ) β . Using these parameters we construct a curve [w 0 /w chiral 0 ](β) and fit it to a quadratic function constrained to pass through 1 for y = 0 There are in total 8 fit parameters, c 1 , c 2 , and six different (a/w chiral 0 ) β . In order to estimate the uncertainty associated with these fit parameters, standard χ 2 fitting does not work since we have statistical uncertainty associated with both w 0 /a and y. The alternative approach that we use is a variant of jackknife resampling. We bin the configurations in steps of 20 to take into account the auto-correlations, which are particularly large for w 0 /a measurements. For each jackknife sample we compute the values of w 0 /a and y. For each sample we perform a χ 2 minimization with χ 2 defined as The error estimate σ (w0/a) is the same for all samples and it is estimated from the fluctuations of this observable over the samples. On the other hand both (w 0 /a) and y change from sample to sample. This minimization produces for each sample eight fit parameters corresponding to c 1 , c 2 , and six values of (a/w chiral 0 ) β . We use these values to estimate fit parameters and their uncertainties. The results are presented in Table II.
The results of the fit are presented in Fig. 2. The measurements for w 0 /a from Table I are rescaled using the (a/w chiral 0 ) β fit parameters from Table II and the solid line represents the chiral extrapolation using the c 1,2 fit results. The dotted line represents the linear chiral extrapolation presented in Ref. [39]. We convert the results in Fig. 4 of Ref. [39],  to straight line in our extrapolation of w 0 as a function of (m π w 0 ) 2 . For this we calculate w 2 0 /t 0 as a function of y; the results are plotted in Fig. 3. For small values of y this ratio can be fit with a straight line and we find the ratio of w 2 0 /t 0 at the chiral point to the value at m π = 390 MeV, the reference point used in Ref. [39], to be 1.31/1.26. Since the ratio of t 0 between these points is 1.08, we find that the ratio of w 0 is 1.06, which leads to a slope of 0.50 since at the reference mass y = 0.11.
Using the functional form w 0 (y)/w chiral 0 = 1−c 1 y+c 2 y 2 we can extrapolate the w 0 /a measurements in Table I to the chiral limit and extract the lattice spacing using the physical value for w chiral 0 . This procedure is used to compute the lattice spacing column in that table and we see that it produces consistent results for the same β but different κ values. Recall that the error in w 0 /a is on the order of a few parts per thousand, so that the overwhelming source of statistical error in the lattice spacing determination comes from w chiral 0 . More accurate determinations of w chiral 0 in physical units will therefore have a significant impact on the accuracy of the lattice spacing determination.

IV. PARAMETERIZATIONS
The results so far are at selected values of β and κ. In this section we perform various interpolations parameterized by smooth functions, to facilitate the usage of the action at other desired lattice spacings and pion masses.

A. Lattice spacing
Using the values for w 0 /a β coming from the global fit described in Section. III, we can perform a quadratic fit to find a parameterization of a as a function of β. The result is log (w 0 /a) = 0.5514 + 1.1600∆β − 0.1217(∆β) 2 , (19) where ∆β ≡ β − 7.3. We use β = 7.3 as the expansion point since the lattice spacing is close to 0.1 fm here. This fit is illustrated in Fig. 4. The error bars are plotted in the figure, but due to the small errors on the a β values, they are not visible.
Using a physical determination for w 0 , this parameterization can be used to obtain the β value corresponding to a desired lattice spacing a in physical units. The function is fairly smooth and we expect that it works reasonably well even outside the range of β used in our simulations.

B. Tadpole factor
The tadpole factor u 0 is an important input parameter that guarantees the improved properties of the Lüscher-Weisz gauge action. Prior to generating our ensembles on the 24 3 × 48 lattice, we did a sensitivity study of pion mass on the input value of u 0 . To this end we generated two 16 3 × 32 ensembles at β = 7.3 and κ = 0.1275, each with 800 configurations (4000 thermalized trajectories separated in steps of five). We used two values for u 0 , 0.868 and 0.870, and found that the pion mass changed from am π = 0.188(2) to 0.179(2), a deviation of 5%. We conclude that changes of u 0 in the fourth significant digit should produce sub-percent level shifts for the pion mass. Both the input value and measured value for u 0 are given in Table I Here w0 is defined in the chiral limit. The parametrization in Eq. 19 describes very well the measurements. The plot above represents the differences between the parametrization and the measured values.
our simulations is very small and we decided to use the same input u 0 for all values of κ.
To interpolate u 0 as a function of β, we use the average measured values of u 0 for the ensembles with the same β to perform a quadratic fit. The result is a smooth function given by u 0 (β) = 0.87010 + 0.03721∆β − 0.01223(∆β) 2 , (20) and displayed in Fig. 5.
This function can be used in conjunction with the function in Eq. 19 to pick the input β and u 0 values corresponding to a desired lattice spacing a.

C. Pion mass
Here we describe the parametrization for the pion mass as a function of the bare parameters β and κ. We performed a series of fits for each set of β values and verified that in the range of quark masses used in our simulations the leading order χPT relation m 2 π ∝ m q describes the data well. We use a linear fit for the fit form to determine the critical value κ c where the pion mass vanishes and the slope. Incidentally the quark mass for Wilson fermions require an additive renormalization; the bare quark mass is then The values of κ c and the slope in the linear fit vary with β. We find that a qudratic fit for the slope and a cubic fit for κ c describe the data well. To determine the parameters, we fit this functional form to all the datapoints in our set. We find the following interpolation function: The interpolation is compared with the data in Fig. 6. We see that the interpolation works well: most of the data points fall within 1% from the curve with the maximal difference at about 3% level. As a check, we used the scale-setting functions to find β, u 0 , and κ values to generate two new ensembles: one on a 28 3 × 56 lattice with a = 0.106 fm and am π = 0.165, and the other one on a 32 3 × 64 lattice with a = 0.093 fm and am π = 0.143 .We measured the lattice spacing on them to be 0.1057(8) fm and 0.0927(7) fm, respectively. The deviation of the central value of our measured lattice spacings from the targeted values are at the sub-percent level. The pion mass am π on these ensembles was measured to be 0.1594(3) and 0.1379(4), with the deviation approximately 4% between the central value of the measured masses and the target values.

V. CONCLUSION
The nHYP smeared, clover improved fermion action with tadpole-improved Lüscher-Weisz gauge fields offers an appealing option for dynamical simulations of hadron physics in QCD. It has good scaling properties, and is relatively inexpensive to simulate when compared to other popular lattice actions such as domain-wall or overlap fermions.
We carried out simulations for 22 different combinations of β and κ on a fixed 24 3 × 48 lattice. For our simulations we ensured the finite-volume corrections are under control (at the sub-percent level on the pion masses considered). We measure the lattice spacing using the Wilson flow parameter w 0 and t 0 and the pion mass on these ensembles. To fix the lattice spacing we use a mass-independent method by matching the value of w 0 in the chiral limit. We obtain smooth interpolating functions for lattice spacing a(β), the tadpole improvement factor u 0 (β), and pion mass dependence on the bare couplings m π (β, κ). Interpolations based on these functions should not deviate by more than a few percent in the parameter region where the lattice spacing is between 0.066 fm and 0.115 fm and pion mass between 207 MeV and 834 MeV. We expect that these interpolations would also work reasonably well outside this range.
As a byproduct of the scale setting method used in this paper, we also map out the behaviour of w 0 as a function of the pion mass. This curve is universal, that is it does not depend on the discretization used to simulate N f = 2 QCD. Our determination agrees well with other calculations [39] and extends it to heavier quark masses.

ACKNOWLEDGMENTS
This work was supported in part by DOE Grant No. DE-FG02-95ER40907. RB is supported in part by the U.S. Department of Energy and ASCR, via a Jefferson Lab subcontract No. JSA-20-C0031. AA gratefully acknowledges the hospitality of the Physics Department at the University of Maryland where part of this work was carried out. The computations were performed on the GWU Colonial One computer cluster and the GWU IMPACT collaboration clusters.