Improved lattice fermion action for heavy quarks
Yong-Gwi Cho, Shoji Hashimoto, Andreas Jüttner, Takashi Kaneko, Marina Marinkovic, Jun-Ichi Noaki, Justus Tobias Tsang
aa r X i v : . [ h e p - l a t ] M a y Prepared for submission to JHEP
UTHEP-673, KEK-CP-321, CERN-PH-TH-2015-074
Improved lattice fermion action for heavy quarks
Yong-Gwi Cho, a Shoji Hashimoto, b,c
Andreas J¨uttner, d Takashi Kaneko, b,c
MarinaMarinkovic, d,e
Jun-Ichi Noaki b and Justus Tobias Tsang d a Graduate School of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571,Japan b High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan c School of High Energy Accelerator Science, The Graduate University for Advanced Studies (Sok-endai), Tsukuba 305-0801, Japan d School of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ,UK e CERN, Physics Department, 1211 Geneva 23, Switzerland
E-mail: [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] Abstract:
We develop an improved lattice action for heavy quarks based on Brillouin-type fermions, that have excellent energy-momentum dispersion relation. The leadingdiscretization errors of O ( a ) and O ( a ) are eliminated at tree-level. We carry out a scalingstudy of this improved Brillouin fermion action on quenched lattices by calculating thecharmonium energy-momentum dispersion relation and hyperfine splitting. We present acomparison to standard Wilson fermions and domain-wall fermions. ontents Charm and bottom quarks have substantially shorter Compton wave-lengths than thetypical length scale of Quantum Chromodynamics (QCD), 1 / Λ QCD . This poses a problemfor numerical simulations of QCD on the lattice. The resolution of the lattice, the latticespacing a , is chosen such that a is sufficiently smaller than 1 / Λ QCD while the entire latticesize L has to be much larger than the length scale of the inverse pion mass, the lightestparticle in the system. For the lattices that can be generated with currently availablecomputational resources, the charm quark mass m c is similar to 1 /a , and the bottomquark mass m b is even larger. This was the motivation for introducing the static or non-relativistic effective theories for heavy quarks, which allow for disentangling the relevantphysical scales in these calculations. The clear scale separation helps in the control ofthe systematics, but the effective theory approaches require an increasing number of extraterms and tuning their associated parameters in order to achieve more precise calculations.As an alternative to the effective heavy quark theories, in this work we perform an extensive– 1 –easibility study of different relativistic approaches to the heavy quark physics from thelattice.Treating heavy quarks on the lattice with the conventional relativistic formulationhas the advantage that the calculation can be made more and more precise as smallerlattice spacings become available. Currently, the finest lattices have 1 /a ≃ m is not much smaller than 1 /a . One successfulexample is the Highly Improved Staggered Quark (HISQ) formulation [1], for which thestaggered fermion formulation is improved by introducing higher dimensional operators,and the leading discretization error is of order ( am ) for heavy quarks. This formulationhas been applied for a number of calculations of phenomenologically important quantities,such as D ( s ) and B ( s ) meson decay constants and other form factors [2–6].Among other relativistic actions, which do not involve the complication due to thefermion doubling of the staggered fermion formulations, the widely used formulations stillcontain the discretization effects of O ( a ), which have to be eliminated to achieve a similarlevel of precision to that of the HISQ formulation. This can be done in a systematic wayaccording to the recipe of the Symanzik improvement program [7, 8], and some attemptswere made in the past [9–12] but they have not been used extensively except for the minimalone, i.e. the O ( a )-improved (or clover) action [9], mainly because the non-perturbativetuning of improvement parameters requires a lot of effort.The goal of this work is to study the scaling of relativistic heavy-quark formulationsin the quenched approximation, before dynamical configurations with similar parametersbecome available. In particular, we present the scaling study of heavy-heavy meson correla-tors, while the scaling of the heavy-strange systems will be presented in a future publication.In this paper, we mainly describe a study of the fermion formulation based on theimproved covariant derivative and Laplacian operators [13]. We compare this Brillouinfermion formulation to more standard lattice fermions, such as the non-improved Wilsonfermion formulation and the M¨obius domain Wall fermions (non-smeared and smeared)[14]. In order to investigate the scaling towards the continuum limit, we generate latticegauge ensembles in the range of 1 /a = 2.0–5.6 GeV in the quenched approximation andperform the measurements of heavy-heavy correlators.Among many options explored in [13], we consider a combination of the “isotropic”covariant derivative (iso) and “Brillouin” Laplacian (bri). This so-called Brillouin fermionis designed such that the violation of four-dimensional rotational symmetry is minimized.By such modification, it turned out that the energy-momentum dispersion relation of amassless fermion is much closer to the continuum one compared to that of a standardWilson fermion [13]. In the context of the Symanzik improvement, this is not obvious sincethe leading discretization error of O ( a ) remains with this prescription. But, as far as thetree-level dispersion relation is concerned, the improvement seems to be achieved includinghigher orders of the lattice spacing a . Once the dispersion relation is improved, one canexpect that interaction terms are also improved, since the form of fermion and gauge– 2 –eld interaction is highly constrained in the gauge theory. Namely, one simply replacesthe tree-level derivative terms by the corresponding covariant derivatives by inserting thegauge links.In this work, we consider a further improvement of the Brillouin-based fermion formu-lation according to the Symanzik improvement program. We design the lattice action suchthat the discretization effects of O ( a ) and O ( a ) are eliminated at the tree-level. With ourchoice we find that the continuum-like energy-momentum dispersion relation is satisfiedvery precisely for quark masses up to am ∼ . O ( a ).The mentioned properties of the Brillouin-type fermions are not guaranteed to besatisfied beyond tree-level, and a non-perturbative study is needed to test the size of thescaling violations in the interacting case. In this work we explicitly check the scalingtowards the continuum limit by taking some basic non-perturbative quantities, such as theheavy-meson dispersion relation and hyperfine splitting.This paper is organized as follows. In Section 2 we review the construction of theBrillouin-type fermion and study its improvement according to the Symanzik improvementprogram. At tree-level, we compare the energy-momentum dispersion relation and complexeigenvalue spectrum of the Dirac operator of various formulations. The improved Brillouinfermion has a limitation on the values of quark mass due to a violation of the reflectionpositivity property as discussed in Section 3. Section 4 describes a non-perturbative scalingstudy of the improved Brillouin fermion and its comparison to the standard Wilson fermionand domain-wall fermions. We then conclude in Section 5. The Brillouin-type covariant derivative and Laplacian operators were introduced in [13].(See also, [18], which introduced similar types of operators in a different context.) Wewrite the lattice Dirac operator as S F = X n,m ψ n D ( n, m ) ψ m , (2.1) D ( n, m ) = X µ γ µ ∇ µ ( n, m ) − a △ ( n, m ) + m δ n,m − c sw X µ<ν σ µν F µν δ n,m , (2.2)– 3 –here ∇ µ ( n, m ) and △ ( n, m ) are the generalized covariant derivative term and Laplacian,respectively. The Sheikholeslami-Wohlert (or clover) term [9] could also be introduced witha coefficient c sw when one introduces the field rotation for the O ( a )-improvement, but wedo not consider this possibility in this paper.For the standard Wilson fermion, the derivative operators are ∇ stdµ ( n, m ) = a ( δ n +ˆ µ,m − δ n − ˆ µ,m ) , (2.3) △ std ( n, m ) = a P µ ( δ n +ˆ µ,m − δ n,m + δ n − ˆ µ,m ) (2.4)at tree-level; the gauge interaction is introduced by promoting the hopping terms δ n ± ˆ µ,m to a covariant derivative including a gauge link. In momentum space, they are given as˜ ∇ stdµ ( p ) = ia sin( p µ a ) = i (cid:18) p µ − a p µ + O ( a ) (cid:19) , (2.5)˜∆ std ( p ) = 2 a X µ (cos( p µ a ) −
1) = − p + O ( a ) . (2.6)The leading discretization effects are ones from ∇ stdµ of O ( a ), as well as those of a △ std ,which is O ( a ). We note that the O ( a ) term of ∇ stdµ violates rotational and Lorentzsymmetry.Among many options proposed in [13], the choice of ∇ isoµ and △ bri leads to the mostcontinuum-like dispersion relation. Their explicit forms are ∇ isoµ ( n, m ) = ρ [ δ n +ˆ µ,m − δ n − ˆ µ,m ] + ρ X ν ( = µ ) [ δ n +ˆ µ +ˆ ν,m − δ n − ˆ µ +ˆ ν,m ]+ ρ X ν = ρ ( = µ ) [ δ n +ˆ µ +ˆ ν +ˆ ρ,m − δ n − ˆ µ +ˆ ν +ˆ ρ,m ]+ ρ X ν = ρ = σ ( = µ ) [ δ n +ˆ µ +ˆ ν +ˆ ρ +ˆ σ,m − δ n − ˆ µ +ˆ ν +ˆ ρ +ˆ σ,m ] , (2.7) △ bri ( n, m ) = λ δ n,m + λ X µ δ n +ˆ µ,m + λ X µ = ν δ n +ˆ µ +ˆ ν,m + λ X µ = ν = ρ δ n +ˆ µ +ˆ ν +ˆ ρ,m + λ X µ = ν = ρ = σ δ n +ˆ µ +ˆ ν +ˆ ρ +ˆ σ,m (2.8)with ( ρ , ρ , ρ , ρ ) = (64 , , ,
1) and ( λ , λ , λ , λ , λ ) = (240 , − , − , − , − µ, ν, ρ, σ = ± , ± , ± , ± i.e. µ = ν = ρ = σ . Underthis restriction, these operators connect neighboring lattice sites m in a 3 hypercube with n in its center. By counting hops along the gauge links, they have up to four hops.In order to make these operators gauge covariant, we have to insert gauge links foreach hop. This should be done so that the rotational symmetry under the cubic group isrespected. We average the shortest possible paths in the taxi-driver distance. For two-hopterms there are two paths; three-hop terms have six paths. The most complicated four-hop– 4 –erms have 24 shortest paths to be averaged. For practical implementation of them, seeAppendix A.In momentum space (at tree level), they have the form˜ ∇ isoµ ( p ) = i a sin( p µ a ) Y ν = µ [cos( p ν a ) + 2]= i p µ (cid:18) −
16 ( p µ a ) (cid:19) Y ν = µ (cid:18) −
12 ( p ν a ) (cid:19) + O ( a ) = ip µ (cid:20) −
16 ( pa ) + O ( a ) (cid:21) (2.9)and ˜ △ bri ( p ) = 4 a "Y µ cos (cid:16) p µ a (cid:17) − = − p + O ( a ) . (2.10)The derivative operator ∇ isoµ has O ( a ) discretization effects which are invariant underrotation, thus the name of “iso”.The Brillouin-type Laplacian (2.8) has the interesting structure that the doublers onthe edges of the Brillouin zone have the same mass. Indeed, for (non-zero) momenta ap µ = ( ± π, ± π, ± π, ± π ), the form in the momentum space (2.10) implies that the inducedmass is always 2 /a . Figure 1 shows ˜ △ std and ˜ △ bri in two-dimensional space. It apparentlyshows that the Brillouin-type Laplacian shows a flat tail at the edge of the Brillouin zone.This can also be seen from the eigenvalue spectrum of the Dirac operator. Figure 2 showsthe eigenvalues of the Dirac operator D ( n, m ) for the Wilson and the Brillouin operatorscalculated on a free gauge field background. The eigenvalues of the Wilson-Dirac operatorplotted on a complex plane show five branches on the real axis, corresponding to thedoublers of masses 0, 2 /a , 4 /a , 6 /a and 8 /a . For the Brillouin operator the doublersare all degenerate at 2 /a . Apart from the real axis, the eigenvalues roughly lie on asingle orbit, very similar to those of the overlap-Dirac operator. It suggests that theoperator is close to the overlap operator and the Ginsparg-Wilson relation is satisfied withgood accuracy at least in the free field case. Among similar lattice fermion formulationswhich involve hopping terms within the 3 hypercube [19–22], the Brillouin fermion isadvantageous for both the continuum-like dispersion relation and the eigenvalue spectrumthat approximately respects the Ginsparg-Wilson relation. One useful measures of the discretization effect is the energy-momentum dispersion relation.It is defined through a pole of the fermion propagator, and takes the form E = p m + p in the continuum theory. For the lattice Dirac operator (2.2), the pole is a solution of (cid:18)
12 ˜ △ ( p ) − ma (cid:19) − X µ (cid:16) ˜ ∇ µ ( p ) (cid:17) = 0 (2.11)– 5 – a ∆ ~ std (p) p a p aa ∆ ~ std (p) -8-7-6-5-4-3-2-1 0 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3-4-3.5-3-2.5-2-1.5-1-0.5 0 a ∆ ~ bri (p) p a p aa ∆ ~ bri (p) -4-3.5-3-2.5-2-1.5-1-0.5 0 Figure 1 . Laplacian operator ˜ △ ( p ) shown in a two-dimensional momentum space ( ap , ap )(Other momentum components are assumed to be zero.) The standard ˜ △ std (left) and Brillouin˜ △ bri (right) are shown. -2-1 0 1 2 0 1 2 3 4 5 6 7 8 WilsonBrillouin
Figure 2 . Eigenvalues of the Dirac operator on a complex plane. They are calculated on a freegauge field background for the Wilson fermion (red circles) and the Brillouin fermion (filled greentriangles). for specific forms of ∇ µ and △ . The poles exist in the Minkowski region that is identifiedby assigning the “energy” E as p = iE . There are more than one poles due to the doublerswhich are heavier than the physical mode by O (1 /a ). In the following we only show thedispersion relation for the physically relevant pole unless otherwise stated.The tree-level dispersion relations are shown in Figure 3 and 4 for Wilson and Brillouinfermions, respectively. As we have an application to heavy fermions in mind, we show theresults for the massive case am = 0 . E = p m + p is shown by a solid line, as well.In the massless limit (left panels), the discretization effect is quite significant for Wilsonfermions beyond | a p | & .
5, while the dispersion relation for the Brillouin fermion closely– 6 – a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 3 . Energy-momentum dispersion relation for Wilson fermions. Two plots show therelation in the massless limit am = 0 (left) and a massive case of am = 0 . | a p | ≡ p ( a p ) in three directions parallel to (1,0,0), (1,1,0) and (1,1,1).Corresponding continuum relation E = p m + p is shown by a solid line. a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 4 . Same as Figure 3, but for the Brillouin fermion. follows that of the continuum theory up to | a p | ≃ .
5. For the massive case, am = 0 . | a p | = 0. If we shift the overall energy such that the dispersionrelation agrees with the continuum one as adopted in the non-relativistic effective theoryapproaches, the deviation would become visible above | a p | ∼ a , we obtain– 7 – aE ) ( a p , am ) = (cid:20) ( am ) − ( am ) + 1112 ( am ) −
56 ( am ) (cid:21) + (cid:20) −
23 ( am ) + 76 ( am ) (cid:21) ( a p ) + (cid:20) −
23 + am (cid:21) X i 56 ( am ) (cid:21) + (cid:20) am ) (cid:21) ( a p ) + h ma i X i 1. Itis understood that the solution of the equation (2.11) becomes complex, which is due tothe lack of reflection positivity. It is potentially dangerous since the Wick rotation to theMinkowski space is not doable in such a situation and one has to assume that the reflectionpositivity is recovered if the continuum limit is taken first. It may have a practical problemthat some instability occurs at relatively low momenta, especially for the massive case, aswe discuss in the following sections. So far, we have shown that Brillouin fermion have some advantageous properties, eventhough it still contains the discretization effect of O ( a ). In the following, we attempt toeliminate these leading discretization errors by modifying the action.Since the O ( a ) error of ∇ isoµ keeps the rotational symmetry, its improvement is rela-tively simple. For instance, we may construct an improved Brillouin action as D imp = X µ γ µ (cid:18) − a △ bri (cid:19) ∇ isoµ (cid:18) − a △ bri (cid:19) + c imp a ( △ bri ) , (2.17)– 9 – a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 6 . Same as Figure 3, but for the improved Brillouin action. -2-1 0 1 2 0 1 2 3 4 5 6 7 8 WilsonBrillouinImpBri Figure 7 . Eigenvalues of the lattice Dirac operators on the free background gauge field. Thepoints show the eigenvalues of Wilson (red circles), Brillouin (filled green triangles) and improvedBrillouin (blue diamonds) fermion operators. where we multiply the Laplacian operator from both sides of ∇ isoµ in order to preserve the γ -hermiticity property. The second term is simply squared with an arbitrary (positive)parameter c imp .This form of the improved action resembles the D34 action, but using ∇ isoµ and △ bri as building blocks the energy-momentum dispersion relation is improved. As shown inFigure 6, the dispersion relation gives a good approximation of the continuum up to | a p | ∼ . 5. The Taylor expansion gives( aE ) ( a p , am ) = (cid:2) ( am ) + 2 c imp ( am ) (cid:3) + ( a p ) , (2.18)which has the leading correction of O ( a ) as expected and it does not contain the possibleterm of O ( a ) that violates the rotational symmetry. This is because the building blocks ∇ isoµ and △ bri themselves reduce the Lorentz violating effects.The eigenvalue spectrum of the Dirac operator on the free background gauge field isshown in Figure 7 for the improved Brillouin fermion (blue) together with those of Wilson– 10 – a E ( a p ) | a p | continuumImpBriImp1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuumImpBriImp1 Figure 8 . Same as Figure 3, but for the improved Brillouin action of reduced numerical cost(2.19), which is shown by blue symbols. The improved Brillouin action of the original form (2.17)is also plotted in green for comparison. (red) and Brillouin (green) fermions. The improved Brillouin eigenvalues form a circlestructure similar to that of the Brillouin operator, but the circle is slightly squashed andpressed on the imaginary axis and approaches the continuum limit where the eigenvaluesare purely imaginary.The improved Brillouin operator D imp defined in (2.17) involves multiple applicationsof △ bri , and therefore is numerically more expensive. Instead, we may consider a lessexpensive operator by using the standard operators for the terms introduced to cancel the O ( a ) errors. Namely, we define D imp = X µ γ µ (cid:18) − a △ std (cid:19) ∇ isoµ (cid:18) − a △ std (cid:19) + c imp a ( △ std ) , (2.19)where △ std is the standard lattice Laplacian operator. The energy-momentum dispersionrelation for this modified operator is shown in Figure 8. Unlike the original improved Bril-louin action (2.17), the departure from the continuum relation is apparent already around a | p | & Since the Brillouin-Dirac operator has an eigenvalue distribution very similar to that ofoverlap fermions as demonstrated in Figure 2, it may be an interesting option to use itas a kernel operator for the overlap-Dirac operator. Projection of eigenvalues to the unitcircle in the complex plane would then require minimal numerical effort, i.e. the order ofthe Chebyshev polynomial or the Zolotarev rational function is relatively lower.Another advantage of the overlap fermion is that the discretization effect of the mass-less Dirac operator is restricted to even powers of a due to its exact chiral symmetry. For– 11 – WilsonBrillouinD imp1 Figure 9 . Same as Figure 7, but for D imp : the improved Brillouin action of reduced numericalcost (2.19) with c imp = 1 / instance, if the standard Wilson-Dirac operator is used as a kernel of the overlap construc-tion, the O ( a ) error of Wilson fermions is eliminated and the leading error becomes O ( a ).If the kernel operator is improved up to O ( a ), then the discretization effect of the corre-sponding overlap operator starts from O ( a ). The massless overlap-Dirac operator can bedefined as D ov (0) = 1 Ra (cid:20) X √ X † X (cid:21) , (2.20)where X is a kernel operator with a large (negative) mass ρ and R is often taken tobe proportional to the unit matrix. Then, D ov satisfies the Ginsparg-Wilson relation { D ov , γ } = RaD ov γ D ov . Introducing a mass, the operator is modified to D ov ( m ) = (cid:18) − am ρ (cid:19) D ov (0) + m. (2.21)It is straight-forward to write down the propagator and solve the pole to obtain theenergy-momentum dispersion relation. With the standard Wilson kernel and ρ = 1, therelation at a p = is ( aE ) = ( am ) + 16 ( am ) . (2.22)This implies that the leading discretization effect is indeed O ( a ). For finite momenta,we plot the dispersion relation in Figure 10. One can see that the dispersion relation isvery similar to that of the kernel operator, which is in this case the Wilson-Dirac operator,shown in Figure 3. With the Brillouin operator as a kernel, the dispersion relation isimproved as shown in Figure 11.Improving the overlap fermion action beyond the O ( a ) discretization effects, one hasto modify the construction of the overlap operator of (2.20), because the Ginsparg-Wilson– 12 – a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 10 . Same as Figure 3, but for the overlap fermion action with the standard Wilson kernelat ρ = 1. a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 11 . Same as Figure 3, but for the overlap fermion action with the Brillouin kernel at ρ = 1. relation of the form { D ov , γ } = RaD ov γ D ov (with a constant R ) already includes O ( a )effects. A possible modification is [23] D impov ( m ) = m + (cid:18) − am ρ (cid:19) D ov + a ρ D † ov D ov , (2.23)where D ov is that of (2.20). In order to eliminate the O ( a ) effects, it has to be usedwith an O ( a ) improved kernel operator. We calculate the dispersion relation for thisimproved overlap-Dirac operator with the improved Brillouin operator as a kernel. Theresult is shown in Figure 12, where we observe that a good approximation for the dispersionrelation is maintained up to | a p | ∼ π/ 2. At zero spatial momentum, the relation E = m issatisfied up to an error of O ( a ).The overlap operator (2.20) is usually constructed using a rational approximation,which is numerically expensive. The number of terms to be included in the rational ap-proximation depends on the range of eigenvalues to be treated and on the desired precision.When the kernel operator is already close to the overlap operator as in the case with theBrillouin operator, it is expected that a minimal order of the rational function wouldachieve a sufficient level of approximation. This property is still to be confirmed with– 13 – a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 a E ( a p ) | a p | continuum(1,1,1)(1,1,0)(1,0,0) Figure 12 . Same as Figure 3, but for the improved overlap fermion action with the improvedBrillouin kernel at ρ = 1. actual numerical calculations. Although the advantage of the Brillouin-type Dirac operators is clear, its numerical cost issubstantially higher than that of the standard (or improved) Wilson fermion action. This issimply because the isotropic derivative and the Brillouin Laplacian involves an interactionto 3 − × γ matrix combination 1 ± γ µ works only for theWilson fermion. Therefore, we expect at least twenty times larger computational costs forthe Brillouin operator, and in practice it is several times more, especially when we use the O ( a )-improved version in (2.17). Therefore, in practical applications the improved Bril-louin fermion could be used only for heavy quarks, for which the fermion matrix inversioncan be carried out with small number of conjugate gradient iterations.In Figure 13 we compare the numerical costs for various Dirac operators by measuringthe elapsed time to solve the heavy quark propagator corresponding to the pseudo-scalarmeson mass m P S ≃ × 32 lattice of 1 /a ≃ D † ( m ) D ( m ).From Figure 13 we can see that the Brillouin fermion takes only five-times more time thanWilson fermions does, despite the above expectation. Likewise, the improved Brillouinfermion is only ten-times slower than Wilson fermions. For the Brillouin fermion, twoimplementations are attempted, i.e. the overall smearing strategy (OSS) and the recursiveformula (Rec) as described in the Appendix A. The OSS implementation has an additionalcost, which we did not account for here, due to an uncounted cost to setup diagonal gaugelinks.We also notice that the performance of the computation on the Blue Gene/Q is differentfor different fermion actions. In our implementation, the number of floating point operationper second (GFlops) per node is about 4.0 for Wilson fermions, while that for other actions– 14 – S e c ond s Cost via CG solver Figure 13 . Numerical cost for various Dirac operators. Elapsed time (in sec) to solve the heavyquark propagator using the conjugate gradient method is plotted. The pseudo-scalar meson mass m P S is roughly tuned to 3.0 GeV. Results are plotted for Wilson fermions, Brillouin fermions (Recand OSS implementations), the improved Brillouin fermion, and the domain-wall fermion. For thedomain-wall fermion, the lattice size in the fifth direction is taken as L s = 8. is around 10, which is compared to the peak performance 200 GFlops, because they aremore compute-intensive. The elapsed time is thus relatively shorter for the actions otherthan Wilson. (In other applications, the JLQCD collaboration uses a highly optimizedcode for the Wilson fermion, which performs much better than 15 GFlops depending onthe condition, but for the comparison in Figure 13 we used a more primitive version of theWilson-Dirac operator in order to make a fair comparison. Optimization of the Brillouinoperators is yet to be done.)This relative speed-up of the Brillouin fermion is explained by the number of conjugategradient iterations to converge. Figure 14 shows the squared norm of the residual vectorat every conjugate gradient iteration steps for Wilson fermions, Brillouin fermions, theimproved Brillouin fermion and the domain-wall fermion. The number of iterations isclearly smaller for the Brillouin-type fermions by more than a factor of two. This explainswhy the Brillouin-type fermions are not as slow as we naively expect. It comes from thefact that the largest eigenvalue of | D (0) | is 2 for the Brillouin operator rather than 8 ofWilson fermions. The condition number of the matrix D ( m ) is thus four-times smaller forthe Brillouin-type fermion. In general, higher derivative terms may give rise to some unphysical poles [12], which aresometimes called “ghost” or “lattice ghost.” If the mass of the ghosts are sufficiently large,no physical effect can be observed, but once they come close to the physical pole, ghostsmay distort the physical solution. Practically, it appears as an oscillatory behavior of theEuclidean correlator. For instance, Figure 15 shows the pseudo-scalar meson correlators– 15 – 50 100 150 200 iteration | r | WilsonBrillouinImp. Bri.Domain-wall Figure 14 . Squared norm of the residual vector at every conjugate gradient iteration steps. Datafor Wilson fermions (solid), Brillouin fermions (dotted), improved Brillouin fermion (thin dashed)and the domain-wall fermion (thick dashed) are shown. l og ( C o rr) t/a ma=0.5ma=0.8ma=1.0ma=1.2ma=1.4ma=1.6ma=1.8ma=2.0ma=2.5ma=3.0 Figure 15 . Heavy pseudo-scalar meson correlators for various heavy quark masses. The statisticalerror is smaller than the thickness of lines. Lines show the correlators of different masses: am =0.5, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5 and 3.0 from top to bottom. calculated with improved Brillouin fermions at various quark masses up to am = 3 . 0. For am & am is large.Therefore there is an upper limit on am to avoid such sickness. One has to be careful,because the problem may be hidden even when the resulting correlation function does notshow the oscillatory behavior. – 16 – a E ( a p = ) am q Figure 16 . Energy of the physical pole (lower curve) and of the ghost (upper curve) as a functionof bare quark mass am . At tree-level, we calculate the energy of the physical pole as well as that correspondingto the lightest doubler. Figure 16 shows those for zero spatial momentum as a function of am . The energy from the physical pole follows the expectation E ( | p | = 0) ≃ m , up to am = 0.6–0.7. On the other hand, the doubler mass slightly decreases for larger am and comesclose to the physical pole near am ≃ am ≃ c imp = 1 / 16 instead of c imp = 1 / O (1 /a ),the upper limit of the heavy quark mass would be 0.5–0.6, according to the tree-levelanalysis shown in Figure 16. The effect of the doublers/ghosts on non-perturbative physicalobservables may appear as a larger scaling violation. Such a symptom will be discussed inthe end of the next section. In this section, we describe a non-perturbative scaling study of the improved Brillouinfermion as well as the standard Wilson and domain-wall fermions towards the continuumlimit. We monitor the energy-momentum dispersion relation and hyperfine splitting ofcharmonium-like heavy-heavy mesons on a set of quenched lattices of inverse lattice spacingbetween 2.0 and 5.6 GeV. – 17 – /a β N sep a − [GeV] a − [GeV] L [fm]16 4.41 100 2.00(07) 2.06(04) 1.579(55)24 4.66 200 2.81(09) 2.89(15) 1.686(52)32 4.89 500 3.80(12) 3.81(09) 1.664(51)48 5.20 40,000 5.64(22) N/A 1.683(64) Table 1 . Quenched lattices used in this study. Temporal lattice size is always T /a = 2 L/a . Thefourth column shows a − determined from the gradient flow. The fifth column is an estimate of thelattice scale from the Sommer scale r = 0.49 fm. We generate a set of SU(3) quenched lattices with the tree-level Symanzik gauge actionat β = 4.41, 4.66, 4.89 and 5.20 as summarized in Table 1. All lattices have a roughlyconstant spatial volume L with L ≃ w = 0.176(2) fm [24]. The lattice spacing determined from the Sommer scale r = 0.49 fm is also listed for three coarser lattices. All ensembles are generated with theheatbath algorithm and the measurement is carried out on gauge configurations separatedby N sep heatbath sweeps, so that the auto-correlation can be safely neglected. The numberof statistical samples is around 100 for each β value except for the finest lattice wherewe have 36 independent gauge configurations. The link smearing procedure is appliedon the gauge configurations before using for the measurement of the heavy-heavy mesoncorrelators (except for that of the “unsmeared” domain-wall fermion, as described below).To be explicit, we employ stout smearing [25] with a parameter α = 0 . c imp = 1 / 8. For comparison, we also employ the standard Wilson fermion and theM¨obius domain-wall fermion. For M¨obius domain-wall fermions, we chose two options:smear or unsmear the gauge links. M¨obius domain-wall fermions are essentially the sameas the conventional domain-wall fermions, but they are designed to achieve much smallerviolation of chiral symmetry at a fixed length L s in the fifth dimension [14]. We chose thescale factor b + c to be 2.0 and L s = 8 with the domain-wall height M = 1.0 for thesmeared domain-wall fermion. The residual breaking of chiral symmetry in this setup isfound to be O (1 MeV) [26].In the following we first describe the measurements with the smeared gauge link.Another set of measurements with unsmeared domain-wall fermion is separately discussedbelow. We tune the quark masses so that pseudo-scalar meson masses become close to ourtarget values m P S = 1.0, 1.5, 2.0, 2.5, 3.0 and 3.5 GeV for all three fermion formulations.The numerical results are interpolated to these target values before comparing the finalresults. Since the length of this interpolation is tiny, we only use a linear function betweennearest two data points.We calculate heavy-heavy meson correlators from four different source points in the– 18 –ime direction. We smear the source with a function e − µ smr r with a parameter µ smr tunedfor each mass to obtain better saturation of the ground state. The gauge configurationsare fixed to Coulomb gauge. The mass and smearing parameter are listed in Table 2.The effective masses for the ground state pseudoscalar and vector mesons are shown inFigure 17. The data corresponding to m P S ≃ β = 4.89 latticeare taken as a typical example. We fit the lattice data with a single exponential function(plus the term representing the contribution from the other temporal direction) in a range[ t min , t max ] = [20,32]. To estimate the systematic error due to the fits, we repeat the fitwith larger t min ’s until t min = 29 and take the variation of their central values as the sizeof systematic error. This similar procedure is applied for other β values and for all masses.The fit results and associated error are also shown in the plots by horizontal lines.In the case of the unsmeared M¨obius domain-wall fermion a slightly different strategyand set of parameters are chosen. Since the unsmeared gauge links are relatively rough,we take a larger value of L s ( L s = 12) with M = 1.6. We work with two different sourcetypes (point and complex Z -wall source [27]). For each source type, we also calculate thequark propagator with a gauge-covariant Gaussian smearing applied on either source orsink (or both). We take many quark masses as described in Appendix B. The results ofthe correlator fits are interpolated to the same reference pseudo-scalar masses in Table 2.This interpolation is carried out with two different ansatzes: a linear interpolation betweenthe nearest two and a quadratic interpolation between the nearest three simulated datapoints. The spread of the central values between the two different approaches gives riseto a systematic error that has been taken into account in the analysis of the continuumextrapolation. In the data for the hyperfine splitting we see a variation of the central valueby up to one sigma when varying the fit-range for the two point functions over a widerange. In order to remain on the conservative side we attach a systematic error of the sizeof the statistical error. The effective speed-of-light c eff can be defined as c ( p ) = E ( p ) − E ( ) p , (4.1)which is unity in the continuum theory and therefore gives a useful measure of the violationof Lorentz symmetry.We calculate the energy with lattice momenta a p of (1,0,0), (1,1,0) and (1,1,1) in unitsof 2 π/L , where L is almost the same for our ensembles. Effective masses for these finitemomentum correlators are shown in Figure 18.As mentioned previously, the systematic error from the mass interpolation needs to betaken into account for the case of the unsmeared domain-wall fermion data. We interpolatethe lattice data with various ans¨atze and take the spread of the central values as thesystematic error. We find that this systematic error is subleading to the statistical error inall cases, but is particularly large for the case of the coarsest ensemble. The reason for thislies in the fact that we had to slightly extrapolate to reach m PS = 3.0 GeV, causing the– 19 – P S [GeV] Dirac op. β = 4 . β = 4 . β = 4 . am µ smr am µ smr am µ smr Wilson 1.038 1.12 0.4946 0.29293.5 Imp. Bri. 0.6675 1.12 0.4517 0.56 0.316 0.4DW 0.728 0.38 0.495 0.3446Wilson 0.69 1.0 0.35 0.23.0 Imp. Bri. 0.55 1.0 0.36 0.5 0.25 0.37DW 0.6 0.7 0.398 0.2785Wilson 0.45 0.2102 0.11052.5 Imp. Bri. 0.416 0.84 0.268 0.45 0.184 0.3DW 0.465 0.303 0.2115Wilson 0.2125 0.1 0.02672.0 Imp. Bri. 0.2808 0.8 0.1705 0.4 0.119 0.25DW 0.3305 0.2154 0.149Wilson 0.0361 − − − − − − − Table 2 . Mass parameters given as inputs for each calculation. m P S is a target heavy-heavy(pseudo-scalar) meson mass. The bare mass parameter am is listed for each fermion formulation:Wilson fermions, improved Brillouin fermions (Imp. Bri.), and smeared domain-wall fermions(DW). For unsmeared domain-wall fermions, see Appendix B. The gauge links are smeared in thesemeasurements as described in the text. µ smr stands for a parameter appearing in the exponentialfunction e − µ smr r to define the source. Since the critical mass is not subtracted, the bare mass forWilson fermions (and for Imp. Bri.) can be negative. systematic error to be larger than on the other ensembles. This systematic error is addedin quadrature to the statistical error before performing the continuum limit.The effective speed-of-light thus calculated is plotted as a function of | p | / (2 π/L ) inFigure 19 for the heavy-heavy pseudo-scalar mesons of mass m P S = 1.5 GeV (left) and3.0 GeV (right). The results on the coarsest lattice (at β = 4.41) show substantial deviationfrom the continuum relation c ( p ) = 1 for Wilson and domain-wall fermions. These twoformulations are close to each other as far as the energy-momentum dispersion relation isconcerned as the tree-level analysis suggests. For the lighter meson of m P S = 1.5 GeV,the deviation is already significant 5–10%; for 3.0 GeV, which is close to charmonium,it is greater than 20%. The data for the improved Brillouin fermion are consistent withunity for both 1.5 and 3.0 GeV mesons. We calculated at four other heavy meson massesas listed in Table 2, and found that the results with the improved Brillouin fermion arealways consistent with unity within the error.We show the scaling of the speed-of-light in Figure 20 against the lattice spacing a . The– 20 – a m e ff t/a Wilson PSV 0.78 0.785 0.79 0.795 0.8 0.805 0.81 0.815 0.82 0.825 5 10 15 20 25 30 35 a m e ff t/a Improved Brillouin PSV 0.78 0.79 0.8 0.81 0.82 0.83 0.84 5 10 15 20 25 30 35 a m e ff t/a Domain-wall PSV Figure 17 . Effective mass of the pseudo-scalar (squares) and vector (triangles) mesons forthe mass parameter corresponding to m P S = 3.0 GeV. Data for Wilson fermions (top), improvedBrillouin fermions (middle) and smeared domain-wall fermions (bottom) at β = 4.89 are plotted.Lines show the fit range and fitted value. data for momentum | p | L/ π = 1 are shown; higher momenta are similar but have largererror bars. As one can see, for the heavy-heavy meson of mass m P S = 3.0 GeV, no deviation– 21 – a m e ff t/a Wilson p L /2 π =1p L /2 π =2p L /2 π =3 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 5 10 15 20 25 30 35 a m e ff t/a Improved Brillouin p L /2 π =1p L /2 π =2p L /2 π =3 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.9 5 10 15 20 25 30 35 a m e ff t/a Domain-wall p L /2 π =1p L /2 π =2p L /2 π =3 Figure 18 . Effective mass of the pseudo-scalar meson with p / (2 π/L ) = 1 (squares), 2 (triangles)and 3 (diamonds). Data corresponding to m P S = 3.0 GeV with Wilson fermions (top), improvedBrillouin fermions (middle) and smeared domain-wall fermions (bottom) at β = 4.89 are plotted.Lines show the fit range and fitted value. from the continuum relation c ( p ) = 1 is found with the improved Brillouin fermion in therange of lattice spacing a . c e ff ( p ) | p | (L/2 π ) WilsonImpBriDomain-wall w/ smearingDomain-wall w/o smearing c e ff ( p ) | p | (L/2 π ) WilsonImpBriDomain-wall w/ smearingDomain-wall w/o smearing Figure 19 . Effective speed-of-light calculated on the coarsest lattice ( β = 4.41) as a functionof momentum squared. Data obtained with Wilson fermions (red squares), improved Brillouinfermions (green circles), smeared domain-wall fermions (blue diamonds), and unsmeared domain-wall fermions (filled magenta pentagons) are shown for pseudo-scalar meson masses at 1.5 GeV(left) and 3.0 GeV (right) c ff ( p ) a[GeV -1 ] m PS =3.00[GeV]WilsonImpBriDomain-wall w/ smearingDomain-wall w/o smearing Figure 20 . Scaling of the speed-of-light against a for the heavy-heavy meson mass m P S =3.0 GeV. The results with | p | L/ π = 1 are shown. The different symbols are those of Wilsonfermions (red squares), improved Brillouin fermions (green circles), smeared domain-wall fermions(blue diamonds), and unsmeared domain-wall fermions (filled magenta pentagons). For the detailson the fit curves, see the text. are similar. A substantial deviation of 20–30% is found on the coarsest lattice, whichdecreases to the level of 5% at a ≃ | p | = 2 π/L . Assuming that the continuum relation c ( p ) = 1 isrecovered in the continuum limit, we employ an ans¨atz f ( a ) = 1 + c a + c a for Wilsonfermions. For smeared and unsmeared domain-wall fermions, an ans¨atz f ( a ) = 1 + c a +– 23 – a is used instead, because the O ( a ) and O ( a ) terms are forbidden by chiral symmetry.Strictly speaking, there might be O ( a ) and O ( a ) discretization effects because of non-zero m res , but we assume that the residual breaking of chiral symmetry is negligible in oursetup. For the improved Brillouin fermion, the leading discretization effects are those of O ( a ). We therefore assume a function f ( a ) = 1 + c a .For each fermion formalism, we obtain a reasonable quality of fit with χ / d.o.f . . c = 0.03(11) GeV and c = − for Wilson fermions, c = − and c = − for smeared domain-wall fermions. For theunsmeared one, we obtain c = − and c = − . Also, c =0.03(9) GeV for improved Brillouin fermions at | p | = 2 π/L is obtained. These fit resultsare plotted in Figure 20. These results suggest that the coefficients have a reasonable sizeof O (1 GeV) or less. The coefficient for the improved Brillouin fermion is essentially zeroeven at the order of a .Since improved Brillouin fermions are designed to achieve O ( a ) and O ( a ) improvementonly in the free theory, there is a possibility that significant contributions of O ( a ) and O ( a ) appear due to radiative corrections. A naive order-counting suggests that their sizeis O ( α s a ) or O ( α s a ), respectively. Assuming α s ∼ The hyperfine splitting m V − m P S is also an interesting quantity to investigate scalingviolations. In the non-relativistic effective theory, it arises from the Pauli term of the form ψ † σ · B ψ . As this term has the same form as the clover term of the O ( a )-improved action,it is expected that the hyperfine splitting is very sensitive to the O ( a ) discretization effectsand also possibly to higher order effects.We show the scaling of m V − m P S in Figure 21 at m P S = 3 . a ≃ O ( a )effects. With domain-wall fermions (both smeared and unsmeared), we also find a signifi-cant discretization effect, though much less severe than in the case of Wilson fermions. Incontrast, the result with improved Brillouin fermions shows very mild a dependence. Fromtheir results alone, one cannot tell any sign of the discretization effects.Here again, we examine the scaling violation by fitting the lattice data as a function of a . Since the value in the continuum limit is unknown, we assume that four lattice formula-tion give the universal result in the continuum limit and fit all the data simultaneously. For– 24 – m V - m PS [ G e V ] a[GeV -1 ] m PS =3.00[GeV]WilsonImpBriDomain-wall w/ smearingDomain-wall w/o smearing Figure 21 . Scaling of the hyperfine splitting against a for the heavy-heavy meson mass m P S =3.0 GeV. The different symbols are those of the Wilson fermion (red squares), the improved Bril-louin fermion (green circles), the smeared domain-wall fermion (blue diamonds) and the unsmeareddomain-wall fermion (filled magenta pentagons). For the details on the fit curves, see the text. m V - m PS [ G e V ] a[GeV -1 ] m PS =3.00[GeV]WilsonImpBriDomain-wall w/ smearingDomain-wall w/o smearing Figure 22 . Same as Figure 21, but with a fit function f ( a ) = 1 + c a + c a for the improvedBrillouin fermion. Wilson fermions we employ f ( a ) = c + c a + c a , while for smeared and unsmeared thedomain-wall fermions we take f ( a ) = c + c a + c a as constrained by chiral symmetry.For improved Brillouin fermions, we first attempt a fit with a function f ( a ) = c + c a . Asshown in Figure 21, a combined fit of four formulations is unsuccessful ( χ /d.o.f. ≃ c = 0.068(1) GeV.It may indicate that the improved Brillouin action receives significant radiative cor-rection at O ( a ) (and O ( a )). Therefore we also try to fit with f ( a ) = c + c a + c a forthe action as that for Wilson fermions. The quality of the fit is better ( χ /d.o.f. ≃ . c is slightly higher: c = 0.077(3) GeV. The results are shown inFigure 22. The fit that allows discretization effects of O ( a ) and O ( a ) indicates that thecoefficients for improved Brillouin fermions ( c = − , c = 0.08(3) GeV ) aremuch smaller than those of Wilson fermions ( c = − , c = 0.18(2) GeV ).One has to be careful about the result for Brillouin fermions, because if this fit captures theactual discretization effect there is a significant cancellation between the O ( a ) and O ( a )terms and the data points between a = 0.25 and 0.5 GeV − show a flat behavior. More de-tailed scaling analysis of other quantities need to be performed if it is the case. For smearedand unsmeared domain-wall fermions, c = − , c = − and c = − , c = − are obtained, respectively. The size of discretizationeffects for these formulations is very similar, though the data point of smeared domain-wallfermion at a = 0.5 GeV − shows significantly larger deviation from the continuum limit.Also we observe that our value of the continuum limit is consistent with that of Ref. [30]in which the charmonium spectra are simulated on quenched lattices with some differentvalence quark formalisms.Finally, we study the scaling of the hyperfine splitting for different sets of heavy quarkmasses, in order to check the discretization effects for the improved Brillouin fermion action.Figure 23 shows the hyperfine splitting for various heavy-heavy meson masses plottedagainst the lattice spacing a . A good scaling is observed in general, but focusing on theheaviest mass ( m P S = 3.5 GeV) we observe a very strong scaling violation for the coarsestlattice point a ≃ − . At this lattice spacing, the natural heavy quark scalingthat the hyperfine splitting decreases for heavier masses is lost between m P S = 3.0 GeVand 3.5 GeV. This may indicate the problem of the improved action for too large am asdiscussed in Section 3. The bare mass for these two masses is 0.55 and 0.67, respectively,which is in the mass range where the effects of doublers could become significant. The energy-momentum dispersion relation is a key property of relativistic field theories.On the lattice, the Lorentz symmetry is explicitly broken by the lattice discretizationand the continuum dispersion relation is expected to be recovered only in the continuumlimit. For practical applications, how fast one can approach the continuum limit becomesa crucial question; while charm quarks can be simulated with moderately large valuesof the input quark mass in lattice units am on ensembles available today, the bottomquark mass cannot be made much less than 0.5–1.0. It is therefore important to design alattice formulation that respects the symmetry relation of the continuum theory as muchas possible. The Brillouin fermion is among such class of fermion formulations, and inthis paper we consider its further improvement and carry out the corresponding numericaltests.The improved Brillouin fermion defined by (2.17) has the energy-momentum dispersionrelation which is close to the continuum one. This is confirmed at the tree-level for themassless case ( am = 0) and the massive case ( am = 0 . O ( a ), but as far as the dispersion relation is concerned, the actual error seems– 26 – m V - m PS [ G e V ] a[GeV -1 ] m PS =1.0 GeVm PS =1.5 GeVm PS =2.0 GeVm PS =2.5 GeVm PS =3.0 GeVm PS =3.5 GeV m V - m PS [ G e V ] a[GeV -1 ] m PS =2.5 GeVm PS =3.0 GeVm PS =3.5 GeV Figure 23 . Hyperfine splitting of the heavy-heavy mesons of mass m P S = 1.0, 1.5, 2.0, 2.5, 3.0,and 3.5 GeV calculated with the improved Brillouin fermion action (top panel). The lattice resultswith the improved Brillouin fermion are plotted as a function of lattice spacing a . The plot in thebottom panel enlarges the results for three heaviest masses. to be much smaller than a naive estimate O (( am ) ) ∼ 13% for am = 0 . 5. Throughthe radiative correction, the terms of O ( aα s ) and O ( a α s ) are induced, and their effectshave to be carefully examined. In our non-perturbative test in the quenched theory, thediscretization effect is insignificant on the lattices of 1 /a = 2.0–5.6 GeV and charm quarkmass m ≃ O ( a ) and O ( a ) effects found in the hyperfinesplitting needs to be studied more carefully, though. The most direct approach would beto calculate on finer lattices. We also to plan to inspect the size of corresponding one-loopcorrection. In successful, the action is a promising candidate for the simulation of charmquark on more realistic unquenched lattices.For more extensive calculations, the numerical cost would become an important issue.With our current implementation, the inversion of charm quark propagators takes 2–3times more time than for the inversion of the domain-wall fermion. It would be worthspending such numerical cost given the highly suppressed discretization effects, but furtherimprovement of the numerical code is certainly desired.– 27 –nother important conclusion from our analysis is that the continuum extrapolationwith the (smeared and unsmeared) domain-wall fermion is possible for charm quark, pro-vided that the data are available in the region beyond 1 /a & References [1] E. Follana et al., Highly improved staggered quarks on the lattice, with applications to charmphysics , Phys.Rev. D75 (2007) 054502, [ hep-lat/0610092 ].[2] E. Follana, C. Davies, G. Lepage, and J. Shigemitsu, High Precision determination of the pi,K, D and D(s) decay constants from lattice QCD , Phys.Rev.Lett. (2008) 062002,[ arXiv:0706.1726 ].[3] C. Davies, C. McNeile, E. Follana, G. Lepage, H. Na, et al., Update: Precision D s decayconstant from full lattice QCD using very fine lattices , Phys.Rev. D82 (2010) 114504,[ arXiv:1008.4018 ].[4] C. McNeile, C. Davies, E. Follana, K. Hornbostel, and G. Lepage, High-Precision f B s andHQET from Relativistic Lattice QCD , Phys.Rev. D85 (2012) 031503, [ arXiv:1110.4510 ].[5] H. Na, C. T. Davies, E. Follana, G. P. Lepage, and J. Shigemitsu, The D → K, lν Semileptonic Decay Scalar Form Factor and | V cs | from Lattice QCD , Phys.Rev. D82 (2010)114506, [ arXiv:1008.4562 ].[6] H. Na, C. T. Davies, E. Follana, J. Koponen, G. P. Lepage, et al., D → π, lν SemileptonicDecays, | V cd | and 2 nd Row Unitarity from Lattice QCD , Phys.Rev. D84 (2011) 114505,[ arXiv:1109.1501 ].[7] K. Symanzik, Continuum Limit and Improved Action in Lattice Theories. 1. Principles and φ Theory , Nucl.Phys. B226 (1983) 187.[8] K. Symanzik, Continuum Limit and Improved Action in Lattice Theories. 2. O(N) NonlinearSigma Model in Perturbation Theory , Nucl.Phys. B226 (1983) 205. – 28 – 9] B. Sheikholeslami and R. Wohlert, Improved Continuum Limit Lattice Action for QCD withWilson Fermions , Nucl.Phys. B259 (1985) 572.[10] T. Eguchi and N. Kawamoto, Improved Lattice Action for Wilson Fermion , Nucl.Phys. B237 (1984) 609.[11] H. W. Hamber and C. M. Wu, Some Predictions for an Improved Fermion Action on theLattice , Phys.Lett. B133 (1983) 351.[12] M. G. Alford, T. Klassen, and G. Lepage, Improving lattice quark actions , Nucl.Phys. B496 (1997) 377–407, [ hep-lat/9611010 ].[13] S. Durr and G. Koutsou, Brillouin improvement for Wilson fermions , Phys.Rev. D83 (2011)114512, [ arXiv:1012.3615 ].[14] R. C. Brower, H. Neff, and K. Orginos, The M¨obius Domain Wall Fermion Algorithm , arXiv:1206.5214 .[15] P. H. Ginsparg and K. G. Wilson, A Remnant of Chiral Symmetry on the Lattice , Phys.Rev. D25 (1982) 2649.[16] H. Neuberger, Exactly massless quarks on the lattice , Phys.Lett. B417 (1998) 141–144,[ hep-lat/9707022 ].[17] H. Neuberger, More about exactly massless quarks on the lattice , Phys.Lett. B427 (1998)353–355, [ hep-lat/9801031 ].[18] M. Creutz, T. Kimura, and T. Misumi, Index Theorem and Overlap Formalism with Naiveand Minimally Doubled Fermions , JHEP (2010) 041, [ arXiv:1011.0761 ].[19] W. Bietenholz and I. Hip, The Scaling of exact and approximate Ginsparg-Wilson fermions , Nucl.Phys. B570 (2000) 423–451, [ hep-lat/9902019 ].[20] W. Bietenholz, Approximate Ginsparg-Wilson fermions for QCD , hep-lat/0007017 .[21] C. Gattringer, A New approach to Ginsparg-Wilson fermions , Phys.Rev. D63 (2001) 114501,[ hep-lat/0003005 ].[22] C. Gattringer, I. Hip, and C. Lang, Approximate Ginsparg-Wilson fermions: A First test , Nucl.Phys. B597 (2001) 451–474, [ hep-lat/0007042 ].[23] H. Ikeda and S. Hashimoto, O( a ) improvement of the overlap-Dirac operator , PoS LAT2009 (2009) 082, [ arXiv:0912.4119 ].[24] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, et al., High-precision scale settingin lattice QCD , JHEP (2012) 010, [ arXiv:1203.4469 ].[25] C. Morningstar and M. J. Peardon, Analytic smearing of SU(3) link variables in latticeQCD , Phys.Rev. D69 (2004) 054501, [ hep-lat/0311018 ].[26] S. Hashimoto, S. Aoki, G. Cossu, H. Fukaya, T. Kaneko, et al., Residual mass infive-dimensional fermion formulations , PoS LATTICE2013 (2014) 431.[27] M. Foster and C. Michael, Quark mass dependence of hadron masses from lattice QCD , Phys.Rev. D59 (1999) 074503, [ hep-lat/9810021 ].[28] A. X. El-Khadra, A. S. Kronfeld, and P. B. Mackenzie, Massive fermions in lattice gaugetheory , Phys.Rev. D55 (1997) 3933–3957, [ hep-lat/9604004 ].[29] T. A. DeGrand, One loop matching coefficients for a variant overlap action: And some of itssimpler relatives , Phys.Rev. D67 (2003) 014507, [ hep-lat/0210028 ]. – 29 – 30] S. Choe et al., Quenched charmonium spectrum , JHEP (2003) 022, [ hep-lat/0307004 ].[31] C. Lehner, Automated lattice perturbation theory and relativistic heavy quarks in theColumbia formulation , PoS LATTICE2012 (2012) 126, [ arXiv:1211.4013 ].[32] G. Cossu, J. Noaki, S. Hashimoto, T. Kaneko, H. Fukaya, et al., JLQCD IroIro++ latticecode on BG/Q , arXiv:1311.0084 . – 30 – Implementation of gauge invariant operators In order to keep the rotational symmetry under the cubic group, one has to average overpossible paths connecting the interacting points. The most economical way is to take thepaths of minimum taxi-driver distance. The average can be represented as off-axis linkvariables. For instance, for the interaction in a ˆ µ -ˆ ν plane connected by two links we maydefine V ˆ µ +ˆ ν ( n ) = 12 [ U µ ( n ) U ν ( n + ˆ µ ) + U ν ( n ) U µ ( n + ˆ ν )] ,V − ˆ µ − ˆ ν ( n ) = 12 h U † µ ( n − ˆ µ ) U † ν ( n − ˆ µ − ˆ ν ) + U † ν ( n − ˆ ν ) U † µ ( n − ˆ µ − ˆ ν ) i ,V ˆ µ − ˆ ν ( n ) = 12 h U µ ( n ) U † ν ( n + ˆ µ − ˆ ν ) + U † ν ( n − ˆ ν ) U µ ( n − ˆ ν ) i ,V − ˆ µ +ˆ ν ( n ) = 12 h U ν ( n ) U † µ ( n + ˆ ν − ˆ µ ) + U † µ ( n − ˆ µ ) U ν ( n − ˆ µ ) i . (A.1)For 3-hop and 4-hop terms, there are off-axis link variables like V ˆ µ +ˆ ν +ˆ ρ ( n ) = 13 [ V ˆ µ +ˆ ν ( n ) U ρ ( n + ˆ µ + ˆ ν ) + V ˆ µ +ˆ ρ ( n ) U ν ( n + ˆ µ + ˆ ρ )+ V ˆ ν +ˆ ρ ( n ) U µ ( n + ˆ ν + ˆ ρ )] ,V ˆ µ +ˆ ν +ˆ ρ +ˆ σ ( n ) = 14 [ V ˆ µ +ˆ ν +ˆ ρ ( n ) U σ ( n + ˆ µ + ˆ ν + ˆ ρ )+ V ˆ µ +ˆ ν +ˆ σ ( n ) U ρ ( n + ˆ µ + ˆ ν + ˆ σ )+ V ˆ µ +ˆ ρ +ˆ σ ( n ) U ν ( n + ˆ µ + ˆ ρ + ˆ σ )+ V ˆ ν +ˆ ρ +ˆ σ ( n ) U µ ( n + ˆ ν + ˆ ρ + ˆ σ )] . (A.2)Note that the Brillouin fermion has 80 nearest-neighbors within a 3 hypercube, and thus80 off-axis link variables have to be prepared. Calculation of these off-axis links can bedone before the conjugate gradient iterations start, as the link variable is unchanged duringthe solver steps. This method is called “overall smearing strategy (OSS)” [13].We also consider an alternative method to calculate the covariant Brillouin Laplacianand the isotropic derivative operators with gauge fields. We use the following recursivedefinition. a △ bri ( n, m ) ψ m = 164 X µ D + µ ψ ′′′ n − ψ n ,ψ ′′′ n ≡ ψ n + 12 X ν = µ D + ν ψ ′′ n ,ψ ′′ n ≡ ψ n + 13 X ρ = µ,ν D + ρ ψ ′ n ,ψ ′ n ≡ ψ n + 14 X σ = µ,ν,ρ D + σ ψ n , (A.3)– 31 –here D ± µ is defined by D ± µ ψ n = U µ ( n ) ψ n +ˆ µ ± U † µ ( n − ˆ µ ) ψ n − ˆ µ . (A.4)We can write down a similar formula for the isotropic derivative. That is ∇ isox ( n, m ) ψ m = 1432 D − x ξ ′′′ n + 12 X ν = x D + ν η ′′′ n ,ξ ′′′ n ≡ ψ n + 12 X ν = x D + ν ξ ′′ n ,ξ ′′ n ≡ ψ n + 13 X ρ = x,ν D + ρ ξ ′ n ,ξ ′ n ≡ ψ n + 14 X σ = x,ν,ρ D + σ ψ n ,η ′′′ n ≡ D − x ξ ′′ n + 13 X ρ = x,ν D + ρ η ′′ n ,η ′′ n ≡ D − x ξ ′ n + 14 X σ = x,ν,ρ D + σ η ′ n ,η ′ n ≡ D − x ψ n . (A.5)This recursive formula gives the same result as OSS does. This compact form is also usefulfor perturbative calculation using the automatic calculation package such as [31].Computational codes for both options are implemented in the Iroiro++ package [32].– 32 – Simulated parameters for the unsmeared M¨obius domain-wall fermion L/a am smearingstart step end parameter16 0.1 0.05 0.4 2.824 0.1 0.05 0.4 4.032 0.1 0.05 0.4 5.648 0.04 0.04 0.28 11.716 0.1 0.05 0.4 4.524 0.066 0.033 0.396 6.032 0.07 0.04 0.39 8.848 0.04 0.04 0.28 11.7 Table 3 . Simulated bare quark mass in lattice units for the unsmeared M¨obius domain-wallfermion. The masses starting from “start” with a step of “step” and ending at “end” are simulated.The “smearing parameter” refers to the choice of the smearing parameter for the Gaussian smear-ing of the source/sink of the propagators. The first block corresponds to the dispersion relationmeasurements with a point source, the second block to the hyperfine splitting measurements with Z -wall source. All measurements we carried out with the unsmeared M¨obius domain-wall fermionare with parameters L s = 12 and M = 1.6.= 1.6.