Scale setting for N f =3+1 QCD
WWUB/20-00
Scale setting for N f = 3 + 1 QCD
LPHA A Collaboration
Roman Höllwieser, Francesco Knechtli, Tomasz Korzec
Dept. of Physics, University of Wuppertal, Gaußstrasse 20, 42119 Germany
Abstract
We present the scale setting for a new set of gauge configurations generated with N f = 3 + 1 Wilson quarks with a non-perturbatively determined clover coefficient in a massive O( a ) improve-ment scheme. The three light quarks are degenerate, with the sum of their masses being equal toits value in nature and the charm quark has its physical mass. We use open boundary conditions intime direction to avoid the problem of topological freezing at small lattice spacings and twisted-mass reweighting for improved stability of the simulations. The decoupling of charm at low energyallows us to set the scale by measuring the value of the low-energy quantity t (cid:63) /a , which is theflow scale t at our mass point, and comparing it to an N f = 2 + 1 result in physical units. Wepresent the details of the algorithmic setup and tuning procedure and give the bare parameters ofensembles with two lattice spacings a = 0 . fm and a = 0 . fm. We discuss finite volumeeffects and lattice artifacts and present physical results for the charmonium spectrum. In particularthe hyperfine splitting between the η c and J/ψ mesons agrees very well with its physical value. a r X i v : . [ h e p - l a t ] A p r ontents Introduction
In [1] a lattice action was proposed to simulate Quantum Chromodynamics (QCD) with N f =3 + 1 quarks on the lattice. The renormalization and improvement conditions are imposed atfinite quark masses. The masses of the up, down and strange quarks are degenerate and theirsum is approximately equal to its value in nature. This not only reduces the number of massparameters, but also facilitates the simulations not having very light up and down quarks. Themass of the charm quark is close to its physical value. Quantities like the charmonium spectrumare fairly insensitive to our approximation of degenerate light quark masses. The action in [1] wasdesigned to simplify the O( a ) ( a is the lattice spacing) improvement pattern for Wilson quarks[2] in the presence of a dynamical charm quark. The improvement coefficients are determinedfollowing the Symanzik improvement program. In a massive scheme for Wilson fermions theSheikholeslami–Wohlert coefficient c sw ( g , aM ) [3] in the action depends, besides on the gaugecoupling g = 6 /β also on the quark masses ( aM being the bare quark mass matrix). The same istrue for the improvement coefficients of operators. The main result of [1] was a non-perturbativedetermination of c sw using a finite volume scheme. In this work we use the action for large volumesimulations to test its scaling properties and provide first physical results for the charmoniumspectrum. The hyperfine splitting m J/ψ − m η c between the ground-state vector and pseudoscalarmesons made from two charm quarks is notoriously sensitive to lattice artifacts stemming fromthe discretization of the Dirac operator [4–6], and therefore a good measure of the quality of theaction.One result of our work is the scale setting for N f = 3 + 1 QCD simulated with the actionof [1]. Scale setting provides a relation between the bare coupling of the simulations and the latticespacing in fm. The scale t (cid:63) (mass dimension -2) is the flow scale t [7] at a particular, unphysicalmass point, namely one where m up = m down = m strange ≡ m l and φ ≡ t (cid:18) m K + m π (cid:19) = 12 t m π, K = 1 . , (1.1) φ ≡ √ t ( m D s + 2 m D ) = √ t m D,D s = 11 . . (1.2)In our setup the pion and kaon as well as D - and D s -mesons are degenerate. In N f = 3 eq. (1.2)plays no role, in our N f = 3 + 1 simulations discussed here, it is important, that the value of φ corresponds to a charm quark mass, about the same as the one found in nature.We measure t ∗ /a at a mass point defined by eq. (1.1) and eq. (1.2) and assume that (cid:112) t (cid:63) =0 . fm to obtain the lattice spacing in fm. This particular value was determined in [8, 9]in N f = 2 + 1 QCD by determining the ratio of t − / , at the unphysical mass point, with alinear combination of pion and kaon decay constants at the physical mass point, in the continuumlimit. Relying on a three flavor result is justified by the fact, that when dealing with low-energyquantities like t , our N f = 3 + 1 theory can be well described by an effective theory whichto leading order is QCD with just N f = 3 light flavors. A detailed study of non-perturbativedecoupling of the charm quark was performed in a model, QCD with N f = 2 mass-degenerateheavy “charm” and zero light quarks [10–13]. The effects of the decoupled charm quark at lowenergies are incorporated in the matching of the gauge coupling and the (less relevant) quarkmasses. Power corrections to decoupling are of the order ( E/m charm ) and (Λ /m charm ) and stemfrom higher dimensional terms in the effective Lagrangian. An important result of the decouplingstudy was that dimensionless low energy quantities are insensitive to the presence of a dynamicalcharm quark at our level of accuracy. Indeed, only effects at the permille level were found, which3s far below the statistical uncertainty of (cid:112) t (cid:63) . Should the precision of (cid:112) t (cid:63) in the N f = 3 theory increase in the future, we might become limited by the accuracy to which decouplingholds. A full-fledged simulation program of N f = 2 + 1 + 1 QCD, including simulations closeto the physical mass point, would then be necessary to reduce the scale-setting error further. Wediscuss the simulation setup and tuning procedure in section 2 and introduce the observables wemeasure in section 3. In section 4 we explain the scale setting of our ensembles with two latticespacings and discuss systematic errors in section 5. These ensembles are not only importantfor a precise determination of the strong coupling constant, but will be very useful for manyother physics applications, e.g. for high precision charm physics. As a first result we thereforepresent charmonium spectra in section 6 and extract charmonium masses and hyperfine splittingin good agreement with PDG [14] values. We end with concluding remarks and a short outlook insection 7.
The ensembles were generated using the open-source (GPL v2) package openQCD version 1.6 [15]in plain C with MPI parallelization. This software has been successfully used in various large-scale projects. The measurements of the hadronic correlation functions were carried out with theopen-source (GPL v2) program “mesons” [16], which is based on various openQCD modules, inparticular the openQCD inverters, and hence shares its architecture (C+MPI) and good scalability.The simulations are performed on lattices of size N t × N s . Gauge fields are periodic in spatialdirections and absent on temporal links sticking out of the lattice, i.e. we have open boundariesin temporal direction imposed on time slice and N t − , hence the physical time extent is T =( N t − a with the lattice spacing a . For a general introduction to lattice QCD simulations withopen boundary conditions and twisted-mass reweighting we refer to [17, 18]. The full action reads S = S g + S f , with a gauge action given by the Lüscher–Weisz action [19,20]with tree-level coefficients c = 5 / and c = − / , S g [ U ] = 1 g (cid:40) c (cid:88) p w ( p )tr [1 − U p ( x )] + c (cid:88) r w ( r )tr [1 − U r ( x )] (cid:41) , (2.1)where the summation is over all oriented plaquettes p and rectangles r (double-plaquettes) on thelattice. U p/r is the product of four/six SU(3) gauge fields U µ ( x ) around the elementary plaquette p or rectangle r . The free parameter of the gauge action is the bare coupling g ≡ /β and theweights w ( p ) = w ( r ) = 1 except for spatial plaquettes and rectangles at T = 0 and T = ( N t − a where w ( p ) = w ( r ) = 1 / , i.e. the coefficient of the gluonic boundary improvement term is setto its tree level value c G = 1 . We do the same with the coefficient of the fermionic boundaryimprovement term c F = 1 . The fermion action then reads S f [ U, ψ, ψ ] = a (cid:88) f =1 (cid:88) x ψ f ( x ) [ D W + m f ] ψ f ( x ) + S SW , (2.2)4ith the Wilson Dirac operator [21] D W = 12 (cid:88) µ =0 { γ µ ( ∇ ∗ µ + ∇ µ ) − a ∇ ∗ µ ∇ µ } , (2.3)where ∇ µ and ∇ ∗ µ are the covariant forward and backward derivatives, respectively, and theSheikholeslami–Wohlert action S SW in eq. (2.2) is needed for O(a) improvement [3].Simulating a physical charm quark results in lattice artifacts due to large cutoff effects oforder proportional to a large value of the lattice charm quark mass ( am c ≈ . . In a massindependent scheme the mass dependent O( a ) lattice artifacts are cancelled by improvement termswhere a mass independent coefficient multiplies aM [22, 23]. Some of these coefficients are onlyknown to one loop in perturbation theory which may be acceptable for the light quarks but not ifthe quark mass is large. Therefore a massive renormalization scheme with close to realistic charmmass m c and Symanzik improvement for Wilson quarks was proposed [1]. This results in mass-dependent renormalization parameters for coupling and quark masses as well as a mass-dependentcoefficient in the clover action S SW = a c sw ( g , aM ) (cid:88) x ¯ ψ ( x ) i σ µν ˆ F µν ( x ) ψ ( x ) , ( ˆ F µν is the standard discretization of the field strength tensor [22]), which we determine using thenon-perturbative fit formula obtained in [1] c sw ( g , aM ) = 1 + Ag + Bg A − . g , A = − . , B = − . . (2.4)While in a full massive scheme the renormalization and improvement factors would dependon all quark masses, in practice a hybrid approach is more feasible, where the dependence onthe heavy quarks (just the charm quark in our case) is absorbed into the factors, but light quarksare treated like in a mass-independent scheme. This allows to use the same value of c sw alongchiral trajectories with constant charm quark mass, while at the same time avoiding large O ( am c ) artifacts. In fact, eq. (2.4) can be used whenever up, down and strange quarks are light, and thecharm quark has close to its physical mass. In this work we describe simulations at the SU (3) flavor symmetric point, where M = diag( m l , m l , m l , m c ) . The generation of these configurations proceeds according to a variant of the Hybrid Monte-Carlo(HMC) algorithm [24]. The classical equations of motion are solved numerically for trajectoriesof length τ = 2 in all simulations, leading to Metropolis proposals which are accepted with anacceptance rate (cid:104) P acc (cid:105) . In order to reduce the computational cost and obtain a high acceptancerate, we split the action and corresponding forces as explained in detail in this section.Since the Wilson Dirac operator is not protected against eigenvalues below the quark mass,we use the second version of twisted mass reweighting, suggested in ref. [25], for the asymmetriceven–odd preconditioned [28] Dirac operator ˆ D = D ee − D eo ( D oo ) − D oe D = D W + m = (cid:18) D ee D eo D oe D oo (cid:19) . (2.5) Often, we quote the hopping parameters κ f = 1 / (2 am f + 8) instead of the bare quark masses. is then simulated with a weight proportional todet [ D † D ] → det [( D oo ) ] det ˆ D † ˆ D + µ ˆ D † ˆ D + 2 µ det [ ˆ D † ˆ D + µ ] (2.6)where we introduced an infrared cutoff by the twisted mass µ . To further reduce the fluctuationsin the forces, we factorize the last term in eq. (2.6) according to [29–31] det[ ˆ D † ˆ D + µ ] = det[ ˆ D † ˆ D + µ N ] × det[ ˆ D † ˆ D + µ ]det[ ˆ D † ˆ D + µ ] × . . . × det[ ˆ D † ˆ D + µ N − ]det[ ˆ D † ˆ D + µ N ] , with aµ i given by { . , . , . , . } for all our ensembles. The individual factors canbe represented by pseudo-fermions such that the combination of twisted-mass reweighting andfactorization leads to a pseudo-fermion action for the light quark doublet with N + 2 = 6 terms S ud , eff [ U, φ , . . . , φ N +1 ] = (cid:0) φ , ˆ D † ˆ D + 2 µ ˆ D † ˆ D + µ φ (cid:1) + N (cid:88) i =1 (cid:0) φ i , ˆ D † ˆ D + µ i ˆ D † ˆ D + µ i − φ i (cid:1) + (cid:8)(cid:0) φ N +1 , D † ˆ D + µ N φ N +1 (cid:1) − D oo } . (2.7)where φ i are six pseudo-fermion fields with support on the even sites of the lattice. The u/d quarkdoublet comes with a reweighting factor W ud = det [ ˆ D † ˆ D ( ˆ D † ˆ D + 2 µ )( ˆ D † ˆ D + µ ) − ] , (2.8)which is estimated stochastically.The strange and charm quarks are simulated with the RHMC algorithm [26,27], decomposingdet ˆ D q = W q det R − q , q ∈ [ s, c ] (2.9)into reweighting factors (to be estimated stochastically again) W q = det [ ˆ D q R q ] (2.10)and Zolotarev optimal rational approximations of ( ˆ D † q ˆ D q ) − / [32] R q = A q N qp (cid:89) k =1 ˆ D † q ˆ D q + ν q,k ˆ D † q ˆ D q + ω q,k (2.11)in the spectral ranges [ r a , r b ] q of ˆ D † q ˆ D q with a number of poles N qp , optimized during the tuningprocess, which uniquely determine the parameters A q , ω q and ν q . The Zolotarev rational func-tion eq. (2.11), in the case of the strange quark, is further broken into several factors, introducingseparate pseudo-fermion fields for the five smallest ω s,k , ν s,k , whereas the determinant of the re-maining factors (first seven poles) is expressed as a single pseudo-fermion integral. In practice,this product is rewritten via partial fraction decomposition into a sum (cid:89) k =1 ˆ D † s ˆ D s + ν s,k ˆ D † s ˆ D s + ω s,k = 1 + (cid:88) k =1 ρ s,k ˆ D † s ˆ D s + ω s,k (2.12) The openQCD code [15] is specialized in simulating a light quark doublet with twisted mass reweighting andfurther quarks with higher masses using the RHMC algorithm [26, 27], that’s why we treat the u/d doublet and strangequark separately, even though they are degenerate. ρ s,k = ( ν s,k − ω s,k ) (cid:89) m =1 ,m (cid:54) = k ν s,m − ω s,k ω s,m − ω s,k . (2.13)The effective action for the strange quark with N sp = 12 poles reads S s , eff [ U, φ , . . . , φ ] = (cid:0) φ , (1 + (cid:88) k =1 ρ s,k ˆ D † s ˆ D s + ω s,k ) φ (cid:1) (2.14) − log det D s oo + (cid:88) l =8 (cid:0) φ l , ˆ D † s ˆ D s + ν s,l ˆ D † s ˆ D s + ω s,l φ l (cid:1) . (2.15)The charm quark contribution to the action is not further factorized S c , eff [ U, φ ] = (cid:0) φ , (1 + N cp (cid:80) k =1 ρ c,k ˆ D † c ˆ D c + ω c,k ) φ (cid:1) − log det D c oo , (2.16)such that we have 13 pseudo-fermion fields and 14 actions in total, with forces summarized intable 9. The molecular dynamics equations are integrated using multi-level higher order (2ndand 4th order OMF) integrators [33]. Linear systems involving the Dirac operators are dealt withusing low-mode-deflation [34] and SAP-preconditioned [35] Krylov space solvers based on localcoherence [36, 37], the tuning of corresponding parameters is discussed in the next section. For the algorithmic parameters, we started with the setup of CLS’s H400 simulation, cf. [38], towhich we added the charm quark. The new contribution to the action was not further factorizedand the corresponding forces were integrated on the intermediate level of our three level integrator.The gauge fields are integrated on the innermost level with the fourth order integrator suggestedby Omelyan, Mryglod, and Folk (OMF) [33], most of the fermion forces are on the intermediatelevel (4th order OMF) and only the most expensive doublet factor (the one with φ in eq. (2.7))and four largest poles of the strange rational function are on the outermost (coarsest) level using asecond order OMF. For most fermion forces we use the locally deflated solver [34, 36, 37] with adeflation subspace block size and 24 deflation modes per block. We set the relative residua to − in the acceptance-rejection step and to − in the force computation. The parameters ofthe rational functions, i.e. number of poles and spectral ranges, are adjusted during thermalizationby calculating the reweighting factors W s,c and the true spectral range of | γ ˆ D | for a subset of thegenerated gauge field configurations. The simulation and algorithmic parameters are summarizedin table 9. The Wilson flow can be used to define quantities with a finite continuum limit in lattice QCD,starting out with a smoothing flow equation [7, 39] ∂ t V µ ( x, t ) = − g { ∂ x,µ S W [ V ] } V µ ( x, t ) , V µ ( x,
0) = U µ ( x ) , (3.1)7ith smeared gauge field V µ ( x, t ) at flow time t and the corresponding Wilson gauge action S W = g (cid:80) p tr [1 − V p ( x )] .Using a symmetrized clover definition of the field strength tensor ˆ G µν ( x, t ) constructed fromthe smooth fields V µ ( x, t ) , the time slice action density E ( x , t ) and the global topological charge Q ( t ) can be defined as E ( x , t ) = − a L (cid:80) (cid:126)x tr { ˆ G µν ( x, t ) ˆ G µν ( x, t ) } , (3.2) Q ( x , t ) = − a π (cid:80) (cid:126)x (cid:15) µναβ tr { ˆ G µν ( x, t ) ˆ G αβ ( x, t ) } . (3.3)Since significant boundary effects in E ( x , t ) and Q ( x , t ) were observed in [38], we take a plateauaverage in a range of x values away from the boundaries (given in table 3) to define the vacuumexpectation value of the energy (cid:104)E ( t ) (cid:105) ant topological charge Q ( t ) . Correlators constructed fromgauge fields at t > are automatically renormalized [40], which allows to define low-energyscales t [7], t c and w [41] as the flow times t at which t (cid:104)E ( t ) (cid:105) = 0 . , . and t (cid:104) ∂ E ( t ) /∂t (cid:105) = 0 . , (3.4)respectively, the partial derivative in the case of w is evaluated numerically. We want to extract meson masses from the zero momentum correlation functions f OO ( x , y ) = a (cid:88) x , y (cid:104)O ( x ) O † ( y ) (cid:105) (3.5)The operators O in table 1 are of the form ¯ q Γ q , where Γ is × spin matrix. Integrating out thefermions leaves us with a single, connected diagram of the form − (cid:88) x , y (cid:10) tr (cid:2) Γ S ( x, y )¯Γ S ( y, x ) (cid:3)(cid:11) gauge . (3.6)with ¯Γ ≡ γ Γ † γ . The trace acts in color and spin space and the propagators S i are given by theinverse Dirac operator (cid:88) y [ D W ( x, y ) + m i ] S i ( y, z ) = δ x,z . (3.7)An efficient way to carry out the calculation is to use stochastic sources, we use 32 noisevectors per time-slice, which amounts to 32 inversions per y value and Dirac structure. In Table 1we list the meson interpolators for the various states investigated and their particle names whichare the closest relatives in nature. An improved signal and exact symmetries are achieved by The notation refers to the γ -structure of the operator. J P C O Particle Plateau A1, A2, BScalar ++ S = ¯ cc (cid:48) χ c {25-43}, {25-43}, {26-43}Pseudoscalar − + P = ¯ uγ d m π {30-70}, {30-90}, {40-100} P = ¯ uγ s m K {30-70}, {30-90}, {40-100} P = ¯ uγ c m D {30-70}, {30-70}, {40-80} P = ¯ sγ c m D s {30-70}, {30-70}, {40-80} P = ¯ cγ c (cid:48) η c {30-70}, {30-70}, {40-80}Vector −− V i = ¯ cγ i c (cid:48) J/ψ {30-70}, {30-80}, {40-100}Axial vector ++ A i = ¯ cγ i γ c (cid:48) χ c {22-35}, {22-35}, {25-35}Tensor + − T ij = ¯ cγ i γ j c (cid:48) h c {18-25}, {18-25}, {22-30} Table 1 : Meson state interpolators and particle names which are the closest relatives in nature.defining the averages ¯ f P P ( x − a ) ≡
12 ( f P P ( x , a ) + f P P ( T − x , T − a )) , (3.8) ¯ f AA ( x − a ) ≡
12 ( f AA ( x , a ) + f AA ( T − x , T − a )) , (3.9) ¯ f V V ( x − a ) ≡ (cid:88) k =1 ( f V k V k ( x , a ) + f V k V k ( T − x , T − a )) , (3.10) ¯ f SS ( x − a ) ≡
12 ( f SS ( x , a ) + f SS ( T − x , T − a )) , (3.11) ¯ f T T ( x − a ) ≡ (cid:88) j>i (cid:0) f T ij T ij ( x , a ) + f T ij T ij ( T − x , T − a ) (cid:1) . (3.12)The meson masses are extracted from the exponential decay of these correlators at x (cid:29) a via weighted plateau averages m = t high (cid:88) x = t low w ( x + a/ m eff ( x + a/ (cid:30) t high (cid:88) x = t low w ( x + a/ , (3.13)of the effective masses am eff ( x + a/ ≡ ln (cid:18) f ( x ) f ( x + a ) (cid:19) . (3.14)and weights w given by their inverse squared errors. Table 1 also lists the plateau ranges for thedifferent ensembles, chosen such that excited state contributions are negligible. With reweighting, QCD expectation values are obtained from the simulated ones via (cid:104)O(cid:105)
QCD = (cid:104)O W (cid:105)(cid:104) W (cid:105) , (3.15)9ith corresponding weights W = W ud W s W c given by the product of reweighting factors eq. (2.8)and eq. (2.10) estimated stochastically. We are interested in functions of primary observables f ( (cid:104)O (cid:105) QCD , . . . , (cid:104)O N (cid:105) QCD , m ) (3.16)and their quark mass derivatives df /dm . Therefore we first reweight our primary observablesaccording to eq. (3.15) and calculate their reweighted mass derivatives via [9] d (cid:104)O i (cid:105) QCD dm = (cid:28) ∂ O i ∂m (cid:29) QCD − (cid:28) O i ∂S∂m (cid:29) QCD + (cid:104)O i (cid:105) QCD (cid:28) ∂S∂m (cid:29)
QCD = (cid:68) ∂ O i ∂m W (cid:69) (cid:104) W (cid:105) − (cid:10) O i ∂S∂m W (cid:11) (cid:104) W (cid:105) + (cid:104)O i W (cid:105) (cid:10) ∂S∂m W (cid:11) (cid:104) W (cid:105) , (3.17)where ∂S∂m q = (cid:80) x ¯ q ( x ) q ( x ) gives rise to (cid:80) x tr [ S ( x, x )] after integration over the quark fields inthe third term of eq. (3.17) and also in the second term, if no other quark fields are present in theobservable O i . This trace is evaluated stochastically using 32 U(1) random sources per configura-tion. When O i contains fermionic fields, the second term in eq. (3.17) produces additional Wickcontractions that depend on O i . For the case of the meson correlators these are worked out andmeasured. For derived observables f , we then use the chain rule d f ( (cid:104)O (cid:105) QCD , . . . , (cid:104)O N (cid:105) QCD , m ) dm = N (cid:88) i =1 ∂f∂ (cid:104)O i (cid:105) QCD d (cid:104)O i (cid:105) QCD dm + ∂f∂m . (3.18)The partial derivatives can be calculated analytically or estimated numerically ∂f∂ (cid:104)O i (cid:105) QCD ≈ (cid:20) f (cid:18) . . . , (cid:104)O i W (cid:105)(cid:104) W (cid:105) + h, . . . (cid:19) − f (cid:18) . . . , (cid:104)O i W (cid:105)(cid:104) W (cid:105) − h, . . . (cid:19)(cid:21) / (2 h ) (3.19)where a suitable h is obtained from the fluctuations of O i , and we do not have explicit derivatives ∂f∂m for the derived observables considered here.In order to correct for mis-tunings of φ eq. (1.1) and φ eq. (1.2), we solve the followingsystem of linear equations to get the mass shifts ∆ m l and ∆ m c for light and charm quarks .
11 = φ + (cid:18) dφ dm u + dφ dm d + dφ dm s (cid:19) ∆ m l + dφ dm c ∆ m c .
94 = φ + (cid:18) dφ dm u + dφ dm d + dφ dm s (cid:19) ∆ m l + dφ dm c ∆ m c and use these to correct arbitrary observables f accordingly f corrected = f + (cid:18) dfdm u + dfdm d + dfdm s (cid:19) ∆ m l + dfdm c ∆ m c . (3.20) For starting parameters we choose a bare coupling β = 3 . , light quark masses given by κ u,d,s =0 . and a charm quark mass by κ c = 0 . , obtained from a rough interpolation of the sim-ulation parameters in [1]. We thermalize on spatially smaller lattices and subsequently double the10ns. Ta × L a β κ l κ c a [fm] Lm (cid:63)π MDUs τ exp A0 × × × × Table 2 : Simulation parameters, lattice sizes and statistics of the three ensembles, the second rowvalues for ensembles A2 and B are shifted to the correct tuning points φ and φ , using the strategydescribed in section 3.3, eq. (3.20).spatial dimensions, using eight copies of the smaller lattices as the new initial configurations. Un-fortunately, the smaller lattices are not helpful for parameter tuning, since we are dealing with largefinite volume effects, as discussed in section 5.2. We measured the flow observables and mesonmasses on a more-or-less thermalized subset of configurations of the final lattice size to evaluate φ and φ . With the starting parameters the tuning point was missed by quite a bit, therefore, wesimulate several other mass points and calculate the derivatives of φ and φ with respect to barequark masses, in order to gradually reach the correct tuning point. This is non-trivial, since themass dependence of t and the meson masses in φ and φ go in opposite directions, see table 10and 11 for results on the ensembles A2 and B. With the final parameters κ u,d,s = 0 . and κ c = 0 . , at β = 3 . we produced two high statistics ensembles A1 and A2 with two differ-ent lattice sizes, and a smaller ensemble A0 on an even smaller volume to study finite size effects.Assuming decoupling of the charm quark, i .e., t (cid:63) | N f =3+1 = t (cid:63) | N f =3 + O (1 /m ) , our valueof t /a ≈ . corresponds to a lattice spacing a ≈ . fm, where we used the N f = 3 result (cid:112) t (cid:63) = 0 . fm from [8, 9]. This makes our smaller volume of ensemble A1 ( L/a = 32 ) L ≈ . fm with m π L = 3 . , which is a bit small, but finite size effects seem to be under control,as the comparison with L/a = 48 shows, see further section 5.2. In a next step we tuned an ensem-ble B at a finer lattice spacing on a × lattice with a bare coupling β = 3 . using the sameprocedure as described above, yielding final parameters κ u,d,s = 0 . and κ c = 0 . , anda lattice spacing a ≈ . fm. This gives a physical lattice extent L ≈ . fm and m π L = 4 . .Simulation parameters are summarized in Table 2, flow measurement and topological chargeresults are presented in Table 3 and corresponding masses and tuning results are shown in Table 4.The second row values for ensembles A2 and B are shifted to the correct tuning points φ and φ ,using the strategy described in section 3.3. The corresponding quark mass shifts a ∆ m l and a ∆ m c are given in Table 5 together with the final lattice spacings. A major challenge in the simulation of QCD with heavy quarks, is the control of lattice artifacts. Inorder to asses their size for our ensembles, and to find out whether they behave like the expectedleading scaling violations, we look at two slightly different definitions of t in eq. (3.4), using11ns. plat. t /a τ int,t t c /a τ int,t c w /a τ int,w Q τ int,Q A0 [22:74] 8.83(23) 10(2) 4.12(9) 9(4) - - 0.83(11) 6(1)A1 [22:74] 7.42(4) 16(7) 3.88(1) 11(4) 10.25(12) 24(12) 1.13(4) 5(1)A2 [22:106] 7.37(2) 27(15) 3.86(1) 17(11) 10.15(6) 25(13) 6.55(18) 3(1)7.43(5) 3.89(2) 10.27(11) 6.77(28)B [30:114] 11.60(6) 35(21) 6.00(2) 32(19) 16.57(17) 29(17) 1.63(7) 8(3)11.62(9) 6.01(3) 16.63(25) 1.61(9)
Table 3 : Flow measurement results and topological charge with integrated autocorrelation times.The second column specifies the plateau range of x values where E ( x , t ) and Q ( x , t ) areaveraged. The second row values for ensembles A2 and B are shifted to the correct tuning points φ and φ , using the strategy described in section 3.3, eq. (3.20).ens. N ms am π,K τ int am D,D s τ int φ φ A0 700 0.310(6) 2.74(92) 0.614(17) 2.24(63) 10.22(90) 15.48(43)A1 1954 0.1141(12) 3.86(82) 0.5232(12) 1.75(39) 1.161(22) 12.098(36)A2 1934 0.1111(4) 4.17(90) 0.5234(4) 1.87(34) 1.092(6) 12.058(17)0.1115(3) 0.5160(16) 1.11 11.94B 2000 0.0896(5) 6.59(81) 0.4135(7) 1.87(17) 1.116(12) 11.950(30)0.0892(4) 0.4128(17) 1.11 11.94
Table 4 : Light meson masses and tuning results, the second row values for ensembles A2 and B areshifted to the correct tuning points φ and φ , using the strategy described in section 3.3.either Wilson’s plaquette discretization of the action density E ( t ) in eq. (3.2), or a more symmetricdefinition based on the symmetric “clover” discretization of the field strength tensor, see [7]. Bothdefinitions of t have to agree in the continuum limit, but may differ at finite lattice spacing,hence, the ratio of t clov0 and t plaq0 has to be one up to cutoff effects. The two values of t arehighly correlated and their ratio can be evaluated very accurately. Figure 1 shows one minusthe ratios for ensembles A2 at β = 3 . and B at β = 3 . (shifted to the correct tuning pointusing eq. (3.20)), supposed to vanish in the continuum limit at a rate proportional to a . Dueto the extremely high precision of this particular observable, the data cannot be described by apure a function shown in the left plot of Figure 1, as that would lead to an unacceptably large χ / dof = 19 . . However, in absolute terms, the deviation from pure a scaling is tiny, only about1.7 permille on the coarser lattice for this particular observable, see the right plot of Figure 1.Given that gradient flow observables in general and the plaquette definition of the action densityin particular are known to have large lattice artifacts (cf. [7]) , we are convinced that for themajority of observables continuum extrapolations linear in a will suffice.12ns. a ∆ m l a ∆ m c a [fm]A2 0.00031(6) -0.0043(9) 0.0536(11)B -0.00001(5) -0.0004(12) 0.0428(7) Table 5 : Quark mass shifts to correct the mistuning of φ and φ and final lattice spacings. Figure 1 : The relative difference of the clover and plaquette definitions of the gradient flow scale t , which vanishes in the continuum limit. The left plot shows a linear fit in a through our dataat two lattice spacings. The right plot shows the tiny deviation from pure a scaling through ourfinest lattice. Next, we study finite volume effects of am π following [42, 43], who propose an analytic scalingformula from chiral perturbation theory ( χ PT) in the “p-expansion" m π ( L ) = m π ( ∞ ) (cid:20) ξ π ˜ g ( Lm π ( ∞ ))2 N f + O ( ξ π ) (cid:21) , ˜ g ( x ) = ∞ (cid:88) n =1 m ( n ) √ nx K ( √ nx ) , ξ π = m π (4 πf π ) . (5.1)We take the pion mass and decay constant in ξ π at the SU(3) flavor symmetrical point determinedon the finest lattice in [9], m π ( ∞ ) our result from ensemble A2, N f = 3 and m ( n ) are integermultiplicities listed in [42] for n = 1 . . . . In the left plot of Fig. 2 we show the analytic χ PTprediction together with our results on ensembles A0, A1 and A2 at the same lattice spacing. Inthe range of pion masses and volumes considered the agreement between the one-loop analyticalprediction and our lattice data is poor, especially for Lm π < . . Finite volume effects seem tobe under control for m π L > , as can be seen by direct comparison of ensembles A1 and A2, seee.g. table 3-8, most values agree within errors.Another source of finite volume effects stems from the finite temporal extent with the openboundaries. Translational invariance is broken and only a portion of the temporal lattice suffi-ciently far away from the boundaries can be used in the averages. The right plot of Fig. 2 showsthese boundary effects for some of the flow observables and the topological charge squared.13 (t ) t c2 E(t c ) t d/dt t E(t)| t=w t E(t ) Figure 2 : Left: finite volume scaling effect of am π represented by the analytic χ PT formula inEq. 5.1 from [42]. Right: Flow measurement plateaus for t E ( x , t ) (black stars), t c E ( x , t c ) (blue crosses), t d/dt t E ( x , t ) | t = w (green circles) and squared topological charge Q ( x , t ) (scaled by t (cid:80) x E ( x , t ) / / (cid:80) x Q ( x , t ) , red plusses) of the lattice with the finest latticespacing (ensemble B). Due to open boundary conditions we average only over time slices x /a =30 − , indicated by the vertical dotted lines for the values given in table 3. How well decoupling of a dynamical charm quark holds has been studied in a sequence of pa-pers [10–13] Based on the outcome of these studies we are optimistic that using the t ∗ value fromthe N f = 3 simulations introduces negligible errors for N f = 3 + 1 QCD with a physical charmquark. The studies were carried out in a model, namely N f = 2 QCD with two degenerate charmquarks and no light quarks. Here we are in the position to investigate decoupling in a more real-istic setting with three degenerate light quarks. We look at dimensionless quantities from ratiosof flow observables (cid:112) t /t c and (cid:112) t /w and attempt continuum limit extrapolations from thecorresponding values from ensembles A2 and B, listed in table 8. In figure 3 we plot the ratiosand their extrapolations and compare with corresponding extrapolations on N f = 3 CLS ensem-bles [38]. The ratios were formed from flow measurements with clover (black) and plaquette (red)definitions of the action density E ( t ) . The continuum extrapolations are performed by constrainedlinear fits which take the correlations into account. The deviations between N f = 3 + 1 and N f = 3 in the continuum are far below 1% and of the order of magnitude found in the modelstudy [11]. We notice that for a reliable continuum extrapolation of the N f = 3 data the com-bination of the two different definitions of the action density and fine lattice spacings (cid:46) . fmare necessary. As we already remarked in the discussion of figure 1 the extremely high precisionof the flow observables exposes cut-off effects which cannot be described by a pure a function.This is particularly visible in figure 3 for flow observables based on the plaquette definition of theaction density. There can be different explanations for the observed behavior of the cut-off effects.One is large a contributions. Another is that while the symmetries of the clover-definition of theaction density allow only for even powers of a in the a -expansion of the observable, the same isnot true for the plaquette definition with open boundary conditions [44]. Finally O ( a ) improvedvariants of both the action density and the flow equations should be employed, when possible [44]and we will investigate these variants in the future.14 N f = 3 N f = 3+1 clover definitionplaquette definition Figure 3 : Comparison of continuum limit extrapolations of (cid:112) t /t c (bottom) and (cid:112) t /w (top) ofour N f = 3 + 1 data (right) with corresponding N f = 3 CLS results (left) from [9, 38] includingthe finest ensemble J500, cf. [45].
When approaching small lattice spacings, the HMC algorithm is known to suffer from criticalslowing down. Autocorrelations grow at least ∝ a − , but have been observed to grow much fasterfor quantities like the topological charge [17]. One method to circumvent the freezing of thetopological charge, is to use open boundary conditions in one of the lattice directions [17], andthat is what we are doing here. As can be seen in Fig. 4 the topological charge moves freely, asdo other “slow” quantities like the flowed action density. Other quantities exhibit the expectedslowing down ∝ a − , which however is still manageable on our finest lattice. E.g. we findintegrated autocorrelation times τ int ,Q ≈ ± [4 MDU] for the topological charge squared and τ int ,t ≈ ± [4 MDU] for t . These values are large, but small enough compared to our totalstatistics. This is corroborated by fits to the tail of the autocorrelation functions which allow usto estimate the exponential autocorrelation time τ exp (corresponding to the slowest mode in theMarkov chain) listed in Table 2. For our error analysis we use the Γ -method [46] adding a tailto the autocorrelation function to account for τ exp [47]. For most observables the errors are onlymarginally different with respect to the standard analysis of [46].15 MDU -10-50510
MDU
YslWsl
Figure 4 : Histories of the topological charge Q ( t ) (left) and of t E ( t ) (right) of the lattice with thefinest lattice spacing (both replica of ensemble B), where t corresponds approximately to t . Thetwo curves show the results for different discretizations of E ( t ) , symmetric “clover” (blue) andplaquette (red), which differ by O ( a ) effects. Due to open boundary conditions we average onlyover time slices x /a = 30 . . . , as indicated in the ordinate label. Although our N f = 3 + 1 ensembles are not at the physical mass point, this work also provides aformidable starting point for N f = 2 + 1 + 1 simulations of QCD with improved Wilson quarks.The quantities φ and φ are chosen such, that the (renormalized) charm quark mass m c andthe sum of the degenerate (renormalized) light quark masses m u + m d + m s are very close totheir physical values. The physical mass point can be approached along chiral trajectories where m u = m d is decreased and m s is increased, while the trace of the quark mass matrix is keptconstant, i.e. the sum of the differences to the physical masses is zero (∆ m u + ∆ m d + ∆ m s = 0) .The chiral extrapolations are expected to be very flat for quantities that do not contain light ( u , d , s )valence quarks. The derivatives of these quantities with respect to the light quark masses areequal to each other at the mass-symmetric point, so in the expression for the correction betweensymmetrical and physical mass point, the leading term vanishes, e.g. for the mass of the η c meson m phys η c = m η c + (∆ m u + ∆ m d + ∆ m s ) dm η c dm l + O (∆ ) , (5.2)because ∆ m u = ∆ m d = − . m s and we only have O(∆ ) corrections. Chiral extrapolationswhere the sum of the light quarks is kept constant have been successfully employed for N f = 2+1 simulations with Wilson fermions [38, 48]. The new ensembles with the fine lattice spacings are very well suited for a study of charmonia,already at the coarsest lattice, the charmonium masses that we measure, neglecting disconnectedcontributions at the moment, are very close to their values in nature. In figure 5 we show themeson spectra of our ensembles A2 and B. We get a clear signal up to the
J/ψ state and can extractreasonable plateau values for higher lying states summarized in table 6. We find good agreementwith PDG [14] data because the states contain only charm valence quarks which in our simulations16ave their physical mass value. Further, we get a very precise result for the charmonium hyperfinesplitting ( m J/ψ − m η ) /m η with perfect agreement to the experimentally known value . , citedby the Particle Data Group, see table 7.We also evaluate a number of other dimensionless quantities from ratios of flow observablesand meson masses, summarized in table 8. For the latter and the charmonium masses in physicalunits we attempt continuum limit extrapolations via linear fits of the corresponding values fromensembles A2 and B. These results are also listed in tables 6-8.ens. η c J/ψ χ c χ c h c am eff A1 0.8175(3) 0.8489(6) 0.943(13) 0.970(19) 0.995(16)A2 0.8180(1) 0.8492(2) 0.9418(94) 0.9744(139) 0.9983(116)0.8114(15) 0.8424((18) 0.935(25) 0.978(36) 0.999(32)B 0.6461(2) 0.6705(3) 0.7508(50) 0.7693(94) 0.793(11)0.6453(20) 0.6699(19) 0.7507(56) 0.768(11) 0.792(12) m eff [GeV] A1 3.010(63) 3.126(66) 3.47(12) 3.57(14) 3.66(14)A2 3.001(62) 3.116(65) 3.46(11) 3.58(12) 3.66(12)2.990(67) 3.104(70) 3.44(16) 3.61(21) 3.68(19)B 2.973(52) 3.086(55) 3.456(79) 3.54(10) 3.65(10)2.973(60) 3.086(62) 3.458(81) 3.54(100) 3.65(11)continuum limit 2.93(18) 3.03(19) 3.45(29) 3.47(33) 3.63(35)2.94(20) 3.05(21) 3.49(36) 3.41(46) 3.60(46)PDG [GeV] 2.9834(5) 3.096900(6) 3.4148(3) 3.51066(7) 3.52538(11) Table 6 : Effective lattice and physical masses of charmonium states η c , J/ψ , χ c , χ c and h c together with continuum limit extrapolations and their PDG [14] values, the second row values forensembles A2 and B are shifted to the correct tuning points φ and φ , for details see section 3.3,eq. (3.20). A2 B cont. ( m J/ Ψ − m η ) /m η Table 7 : Charmonium hyperfine splitting ( m J/ Ψ − m η ) /m η = 0 . (PDG) on various ensembles.The last column gives the continuum limit extrapolations of original and shifted results. We presented the scale setting and tuning of N f = 3 + 1 QCD using a massive renormalizationscheme with a non-perturbatively determined clover coefficient from [1]. We produced two en-sembles A1 and A2 with lattice sizes × and × at a lattice spacing a = 0 . fm,see also [49], and an ensemble B on a × lattice at a finer lattice spacing a = 0 . fm.As a first physics result, we measure the masses of the charmonium states η c , J/ψ , χ c , χ c and h c , which agree with their PDG values. In particular, we can reproduce the charmonium17 igure 5 : Effective masses of the pion/kaon, D - and D s -meson, charmonium states η c , J/ Ψ , χ , χ and h c (from bottom to top) on ensembles A2 (top) and B (bottom).18ns. (cid:112) t /t c (cid:112) t /w m D m π m η m π m J/ Ψ m π m χ c m π m χ c m π m h c m π A1 1.3836(16) 0.851(3) 4.58(4) 7.16(7) 7.44(7) 8.26(13) 8.49(18) 8.71(16)A2 1.3820(8) 0.852(1) 4.71(1) 7.36(2) 7.64(2) 8.48(9) 8.77(13) 8.99(11)1.3827(9) 0.851(2) 4.6267 7.27(1) 7.55(1) 8.38(23) 8.77(33) 8.96(29)B 1.3898(15) 0.837(2) 4.62(2) 7.21(4) 7.49(4) 8.38(6) 8.59(9) 8.86(11)1.3900(21) 0.837(4) 4.6267 7.23(1) 7.51(1) 8.41(6) 8.61(10) 8.88(11)cont.1.4043(37) 0.809(7) 4.45(6) 6.95(10) 7.21(11) 8.22(23) 8.27(33) 8.63(35)1.4043(48) 0.809(9) 4.6267 7.16(3) 7.43(4) 8.47(43) 8.32(64) 8.74(60)
Table 8 : Dimensionless quantities on various ensembles, the second row values for ensembles A2and B are shifted to the correct tuning points φ and φ , using the strategy described in section 3.3,eq. (3.20). Note, there is no error for the shifted m D /m π ratio, since it is given by the exact ratio of φ / √ φ . The last two lines correspond to continuum limit extrapolations of original and shiftedvalues via linear fits of corresponding values from ensembles A2 and B.hyperfine splitting ( m J/ Ψ − m η ) /m η within permille level precision of the experimentally knownvalue . . First continuum limit extrapolations of the measured quantities are attempted.The next steps are to double the statistics of ensemble B and produce a third ensemble C ona × lattice at an even finer lattice spacing a = 0 . fm, in order to make more reliablecontinuum limit extrapolations. For the future, we plan to measure disconnected contributions tothe charmonium masses and extract excited states. A further development of this project is thesimulation of N f = 2 + 1 + 1 QCD close to the physical light quark masses. Another longtermgoal is the determination of the strong coupling α S or equivalently the Λ -parameter in the case offour flavors. The currently most precise result comes from lattice simulations [8], see also [50].One of the remaining uncertainty comes from the use of perturbative decoupling for the matchingof the gauge couplings of QCD with four and three flavors. In [12] this uncertainty was estimatedto be below 1.5% for the ratios of the Λ parameters which is still below the 3.5% precision of the Λ -parameter of [8]. But new techniques like [51] promise a larger precision for α S . The runningof α S in four flavor QCD has been computed in [52]. Such a result, can be combined with thescale setting for N f = 3 + 1 QCD presented in this work to obtain α S fully non-perturbatively infour flavor QCD. Acknowledgements
We thank Rainer Sommer and Ulli Wolff for precious advice on the manuscript and AndreiAlexandru, Stefan Dürr and Stefan Schaefer for valuable discussions. We are grateful to ourCLS colleagues for sharing data used in Fig. 3. We gratefully acknowledge the Gauss Center forSupercomputing e.V. ( ) for funding this project by providing com-puting time on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) underGCS/NIC project ID HWU35. R.H. was supported by the Deutsche Forschungsgemeinschaft inthe SFB/TRR55. 19
Simulation parameters ensemble A1 A2 BID nf31H100 nf31N200 nf31I300Force 0, lvl gauge, 0 gauge, 0 gauge, 0Force 1, lvl TM1-EO-SDET, 1 TM1-EO-SDET, 1 TM1-EO-SDET, 1Force 2, lvl TM2-EO, 1 TM2-EO, 1 TM2-EO, 1Force 3, lvl TM2-EO, 1 TM2-EO, 1 TM2-EO, 1Force 4, lvl TM2-EO, 1 TM2-EO, 1 TM2-EO, 1Force 5, lvl TM2-EO, 2 TM2-EO, 2 TM2-EO, 2Force 6, lvl s-RAT-SDET, 1 s-RAT-SDET, 1 s-RAT-SDET, 1Force 7, lvl s-RAT, 1 s-RAT, 1 s-RAT, 1Force 8, lvl s-RAT, 1 s-RAT, 1 s-RAT, 1Force 9, lvl s-RAT, 2 s-RAT, 2 s-RAT, 2Force 10, lvl s-RAT, 2 s-RAT, 2 s-RAT, 2Force 11, lvl s-RAT, 2 s-RAT, 2 s-RAT, 2Force 12, lvl s-RAT, 2 s-RAT, 2 s-RAT, 2Force 13, lvl c-RAT-SDET, 1 c-RAT-SDET, 1 c-RAT-SDET, 1Level 0,nstep OMF4, 2 OMF4, 2 OMF4, 2Level 1,nstep OMF4, 1 OMF4, 1 OMF4, 1Level 2,nstep OMF2, 8 OMF2, 8 OMF2, 8 κ uds κ c c sw aµ aµ aµ aµ N sp , [ r a , r b ] s
12, [0.001,9.0] 12, [0.001,9.0] 12, [0.002,8.0] N cp , [ r a , r b ] c
10, [0.2,8.0] 10, [0.2,8.0] 8, [0.2,8.0] N traj (cid:104) P acc (cid:105) Table 9 : Parameters of the algorithm: We give the forces used for the gauge and fermion fields withtheir integration levels, the integrators for the different levels and number of steps per trajectoryresp. outer level, the hopping, c sw and twisted-mass parameters, the number of poles N p and theranges used in the RHMC for strange and charm quarks, as well as the total length of the Markovchain and the acceptance rates. 20 Mass derivatives O d O /dm u d O /dm d d O /dm s d O /dm c am π am K am D am D s am η c am J/ Ψ am χ c am χ c am h c t /a t c /a w /a Q φ φ Table 10 : Derivatives with respect to quark masses of pion/kaon, D - and D s -mesons and charmo-nium states η c , J/ψ , χ c , χ c and h c effective masses and flow observables on ensemble A2. O d O /dm u d O /dm d d O /dm s d O /dm c am π am K am D am D s am η c am J/ Ψ am χ c am χ c am h c t /a t c /a w /a Q φ φ Table 11 : Derivatives with respect to quark masses of pion/kaon, D - and D s -mesons and charmo-nium states η c , J/ψ , χ c , χ c and h c effective masses and flow observables on ensemble B.21 eferences [1] ALPHA
Collaboration, P. Fritzsch, R. Sommer, F. Stollenwerk, and U. Wolff, “SymanzikImprovement with Dynamical Charm: A 3+1 Scheme for Wilson Quarks,”
JHEP (2018)025, arXiv:1805.01661 [hep-lat] .[2] K. G. Wilson, “Quarks and Strings on a Lattice,” in New Phenomena in SubnuclearPhysics: Proceedings, International School of Subnuclear Physics, Erice, Sicily, Jul 11-Aug1 1975. Part A , p. 99. 1975. [,0069(1975)].[3] B. Sheikholeslami and R. Wohlert, “Improved Continuum Limit Lattice Action for QCDwith Wilson Fermions,”
Nucl. Phys.
B259 (1985) 572.[4]
QCD-TARO
Collaboration, S. Choe, P. de Forcrand, M. Garcia Perez, Y. Liu,A. Nakamura, I. O. Stamatescu, T. Takaishi, and T. Umeda, “Quenched charmoniumspectrum,”
JHEP (2003) 022, arXiv:hep-lat/0307004 [hep-lat] .[5] Hadron Spectrum
Collaboration, L. Liu, G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas,P. Vilaseca, J. J. Dudek, R. G. Edwards, B. Joo, and D. G. Richards, “Excited and exoticcharmonium spectroscopy from lattice QCD,”
JHEP (2012) 126, arXiv:1204.5425[hep-ph] .[6] Y.-G. Cho, S. Hashimoto, A. Jüttner, T. Kaneko, M. Marinkovi´c, J.-I. Noaki, and J. T.Tsang, “Improved lattice fermion action for heavy quarks,” JHEP (2015) 072, arXiv:1504.01630 [hep-lat] .[7] M. Lüscher, “Properties and uses of the Wilson flow in lattice QCD,” JHEP (2010) 071, arXiv:1006.4518 [hep-lat] . [Erratum: JHEP03,092(2014)].[8] ALPHA
Collaboration, M. Bruno, M. Dalla Brida, P. Fritzsch, T. Korzec, A. Ramos,S. Schaefer, H. Simma, S. Sint, and R. Sommer, “QCD Coupling from a NonperturbativeDetermination of the Three-Flavor Λ Parameter,”
Phys. Rev. Lett. no. 10, (2017)102001, arXiv:1706.03821 [hep-lat] .[9] M. Bruno, T. Korzec, and S. Schaefer, “Setting the scale for the CLS flavorensembles,”
Phys. Rev.
D95 no. 7, (2017) 074504, arXiv:1608.08900 [hep-lat] .[10]
ALPHA
Collaboration, M. Bruno, J. Finkenrath, F. Knechtli, B. Leder, and R. Sommer,“Effects of Heavy Sea Quarks at Low Energies,”
Phys. Rev. Lett. no. 10, (2015)102001, arXiv:1410.8374 [hep-lat] .[11]
ALPHA
Collaboration, F. Knechtli, T. Korzec, B. Leder, and G. Moir, “Power correctionsfrom decoupling of the charm quark,”
Phys. Lett.
B774 (2017) 649–655, arXiv:1706.04982 [hep-lat] .[12]
ALPHA
Collaboration, A. Athenodorou, J. Finkenrath, F. Knechtli, T. Korzec, B. Leder,M. Krsti´c Marinkovi´c, and R. Sommer, “How perturbative are heavy sea quarks?,”
Nucl.Phys.
B943 (2019) 114612, arXiv:1809.03383 [hep-lat] .[13] S. Cali, F. Knechtli, and T. Korzec, “How much do charm sea quarks affect the charmoniumspectrum?,”
Eur. Phys. J.
C79 no. 7, (2019) 607, arXiv:1905.12971 [hep-lat] .2214]
Particle Data Group
Collaboration, M. Tanabashi et al. , “Review of Particle Physics,”
Phys. Rev.
D98 no. 3, (2018) 030001.[15] M. Lüscher, S. Schaefer, “openQCD: Simulation programs for lattice QCD,” 2016. http://luscher.web.cern.ch/luscher/openQCD .[16] T. Korzec, “Calculate meson correlators with openQCD,” 2014. https://github.com/to-ko/mesons .[17] M. Lüscher and S. Schaefer, “Lattice QCD without topology barriers,”
JHEP (2011)036, arXiv:1105.4749 [hep-lat] .[18] M. Lüscher and S. Schaefer, “Lattice QCD with open boundary conditions andtwisted-mass reweighting,” Comput. Phys. Commun. (2013) 519–528, arXiv:1206.2809 [hep-lat] .[19] M. Lüscher and P. Weisz, “On-Shell Improved Lattice Gauge Theories,”
Commun.Math.Phys. (1985) 59.[20] M. Lüscher and P. Weisz, “Computation of the Action for On-Shell Improved Lattice GaugeTheories at Weak Coupling,” Phys. Lett.
B158 (1985) 250.[21] K. G. Wilson, “Confinement of Quarks,”
Phys. Rev.
D10 (1974) 2445–2459. [,319(1974)].[22] M. Lüscher, S. Sint, R. Sommer, and P. Weisz, “Chiral symmetry and O(a) improvement inlattice QCD,”
Nucl. Phys.
B478 (1996) 365–400, arXiv:hep-lat/9605038[hep-lat] .[23] T. Bhattacharya, R. Gupta, W. Lee, S. R. Sharpe, and J. M. S. Wu, “Improved bilinears inlattice QCD with non-degenerate quarks,”
Phys. Rev.
D73 (2006) 034504, arXiv:hep-lat/0511014 [hep-lat] .[24] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, “Hybrid Monte Carlo,”
Phys.Lett.
B195 (1987) 216–222.[25] M. Lüscher and F. Palombi, “Fluctuations and reweighting of the quark determinant onlarge lattices,”
PoS
LATTICE2008 (2008) 049, arXiv:0810.0946 [hep-lat] .[26] A. D. Kennedy, I. Horvath, and S. Sint, “A New exact method for dynamical fermioncomputations with nonlocal actions,”
Nucl. Phys. Proc. Suppl. (1999) 834–836, arXiv:hep-lat/9809092 [hep-lat] . [,834(1998)].[27] M. A. Clark and A. D. Kennedy, “Accelerating dynamical fermion computations using therational hybrid Monte Carlo (RHMC) algorithm with multiple pseudofermion fields,” Phys.Rev. Lett. (2007) 051601, arXiv:hep-lat/0608015 [hep-lat] .[28] T. A. DeGrand, “A Conditioning Technique for Matrix Inversion for Wilson Fermions,” Comput. Phys. Commun. (1988) 161–164.[29] M. Hasenbusch, “Speeding up the hybrid Monte Carlo algorithm for dynamical fermions,” Phys. Lett.
B519 (2001) 177–182, arXiv:hep-lat/0107019 [hep-lat] .2330] M. Hasenbusch and K. Jansen, “Speeding up lattice QCD simulations with clover improvedWilson fermions,”
Nucl. Phys.
B659 (2003) 299–320, arXiv:hep-lat/0211042[hep-lat] .[31] S. Schaefer, “Status and challenges of simulations with dynamical fermions,”
PoS
LATTICE2012 (2012) 001, arXiv:1211.5069 [hep-lat] .[32] N. I. Achiezer,
Theory of approximation . Dover Publications, New York, 1992.[33] I. Omelyan, I. Mryglod, and R. Folk, “Symplectic analytically integrable decompositionalgorithms: classification, derivation, and application to molecular dynamics, quantum andcelestial mechanics simulations,”
Computer Physics Communications no. 3, (2003) 272.[34] M. Lüscher, “Local coherence and deflation of the low quark modes in lattice QCD,”
JHEP (2007) 081, arXiv:0706.2298 [hep-lat] .[35] M. Lüscher, “Solution of the Dirac equation in lattice QCD using a domain decompositionmethod,” Comput. Phys. Commun. (2004) 209–220, arXiv:hep-lat/0310048[hep-lat] .[36] M. Lüscher, “Deflation acceleration of lattice QCD simulations,”
JHEP (2007) 011, arXiv:0710.5417 [hep-lat] .[37] A. Frommer, K. Kahl, S. Krieg, B. Leder, and M. Rottmann, “Adaptive Aggregation BasedDomain Decomposition Multigrid for the Lattice Wilson Dirac Operator,” SIAM J. Sci.Comput. (2014) A1581–A1608, arXiv:1303.1377 [hep-lat] .[38] M. Bruno et al. , “Simulation of QCD with N f = + JHEP (2015) 043, arXiv:1411.3982 [hep-lat] .[39] R. Narayanan and H. Neuberger, “Infinite N phase transitions in continuum Wilson loopoperators,” JHEP (2006) 064, arXiv:hep-th/0601210 [hep-th] .[40] M. Lüscher and P. Weisz, “Perturbative analysis of the gradient flow in non-abelian gaugetheories,” JHEP (2011) 051, arXiv:1101.0963 [hep-th] .[41] S. Borsanyi et al. , “High-precision scale setting in lattice QCD,” JHEP (2012) 010, arXiv:1203.4469 [hep-lat] .[42] G. Colangelo, S. Dürr, and C. Haefeli, “Finite volume effects for meson masses and decayconstants,” Nucl. Phys.
B721 (2005) 136–174, arXiv:hep-lat/0503014[hep-lat] .[43] J. Gasser and H. Leutwyler, “Light Quarks at Low Temperatures,”
Phys. Lett.
B184 (1987)83–88.[44] A. Ramos and S. Sint, “Symanzik improvement of the gradient flow in lattice gaugetheories,”
Eur. Phys. J.
C76 no. 1, (2016) 15, arXiv:1508.05552 [hep-lat] .[45] G. S. Bali, L. Barca, S. Collins, M. Gruber, M. Löffler, A. Schäfer, W. Söldner, P. Wein,S. Weishäupl, and T. Wurm, “Nucleon axial structure from lattice QCD,” arXiv:1911.13150 [hep-lat] . 2446]
ALPHA
Collaboration, U. Wolff, “Monte Carlo errors with less errors,”
Comput. Phys.Commun. (2004) 143–153, arXiv:hep-lat/0306017 [hep-lat] . [Erratum:Comput. Phys. Commun.176,383(2007)].[47]
ALPHA
Collaboration, S. Schaefer, R. Sommer, and F. Virotta, “Critical slowing down anderror analysis in lattice QCD simulations,”
Nucl. Phys.
B845 (2011) 93–119, arXiv:1009.5228 [hep-lat] .[48] W. Bietenholz et al. , “Flavour blindness and patterns of flavour symmetry breaking in latticesimulations of up, down and strange quarks,”
Phys. Rev.
D84 (2011) 054509, arXiv:1102.5300 [hep-lat] .[49] R. Höllwieser, F. Knechtli and T. Korzec, “Scale setting for QCD with Nf=3+1 dynamicalquarks,” arXiv:1907.04309 [hep-lat] .[50]
Flavour Lattice Averaging Group
Collaboration, S. Aoki et al. , “FLAG Review 2019,” arXiv:1902.08191 [hep-lat] .[51]
ALPHA
Collaboration, M. Dalla Brida, R. Höllwieser, F. Knechtli, T. Korzec, A. Ramos,and R. Sommer, “Non-perturbative renormalization by decoupling,” arXiv:1912.06001 [hep-lat] .[52]
ALPHA
Collaboration, F. Tekin, R. Sommer, and U. Wolff, “The Running coupling ofQCD with four flavors,”
Nucl. Phys.
B840 (2010) 114–128, arXiv:1006.0672[hep-lat]arXiv:1006.0672[hep-lat]