Beauty mesons in N f =2+1+1+1 lattice QCD with exact chiral symmetry
NNTUTH-20-505B
Beauty mesons in N f = 2 + 1 + 1 + 1 lattice QCDwith exact chiral symmetry Ting-Wai Chiu
1, 2, 3 Physics Department, National Taiwan Normal University, Taipei, Taiwan 11677, R.O.C. Institute of Physics, Academia Sinica, Taipei, Taiwan 11529, R.O.C. Physics Department, National Taiwan University, Taipei, Taiwan 10617, R.O.C.
Abstract
We present the first study of N f = 2 + 1 + 1 + 1 lattice QCD with domain-wall quarks. The( b , c , s ) quarks are physical, while the ( u , d ) quarks are heavier than their physical masses, withthe pion mass ∼
700 MeV. The gauge ensemble is generated by hybrid Monte Carlo simulationwith the Wilson gauge action for the gluons, and the optimal domain-wall fermion action forthe quarks. Using point-to-point quark propagators, we measure the time-correlation functionsof quark-antiquark meson interpolators with quark contents ¯bb , ¯bc , ¯bs , and ¯cc , and obtain themasses of the low-lying mesons. They are in good agreement with the experimental values, plussome predictions which have not been observed in experiments. Moreover, we also determine themasses of ( b , c , s ) quarks. a r X i v : . [ h e p - l a t ] A p r . INTRODUCTION In 2007, we performed the first study of treating valence ( u , d , s , c , b ) quarks as Diracfermions in quenched lattice QCD with exact chiral symmetry [1, 2]. The low-lying massspectra of mesons with quark contents ¯bb , ¯bc , ¯bs , and ¯cc were determined, together withthe pseudoscalar decay constants. Some of our results (e.g., the masses of η b and h b ) weretheoretical predictions at the time of publication, which turn out to be in good agreementwith later experimental results. This asserts that it is feasible to treat ( u , d , s , c , b ) valencequarks as Dirac fermions, in lattice QCD with exact chiral symmetry.Now the question is whether one can simulate dynamical ( u , d , s , c , b ) quarks in latticeQCD with exact chiral symmetry. This motivates the present study. Since the b quark isheavy, with mass m b ∼ c , it requires a fine lattice spacing such that the condition m b a < L has to be sufficiently large such that M π L (cid:29)
1. These two constraints ( a ∼ .
033 fm and M π L ∼ −
6) together give the lattice size ∼ − (see Fig. 1), which is beyondthe capability of the present generation of supercomputers. Nevertheless, even before the N x M G e V ] L = 2M L = 4M L = 6 M phys b a ~ 0.8 FIG. 1: The design of lattice QCD with physical ( u , d , s , c , b ) quarks. next generation of Exaflop supercomputers will be available ∼ b , c , s ) quarks with physical masses can be dynamically2imulated on the lattice, while keeping u and d quarks heavier than their physical masses.If the pion mass is kept at ∼
700 MeV/ c , then both constraints M π L > m b a < ×
64 lattice. For domain-wall fermion with the extent N s = 16 in thefifth dimension, the entire hybrid Monte Carlo (HMC) simulation [3] on the 40 × × ¯bb , ¯bc , ¯bs and ¯cc . Insection 4, we determine the masses of ( b , c , s ) quarks. In section 5, we conclude with someremarks. II. SIMULATION OF LATTICE QCD WITH EXACT CHIRAL SYMMETRYA. Preliminaries
Since all quarks in QCD are excitations of Dirac fermion fields, it is vital to preserve thisessential feature in lattice QCD. The most theoretically appealing lattice fermion scheme isthe domain-wall/overlap fermion [5–7], which preserves the exact chiral symmetry at finitelattice spacing, thus provides a proper formulation of QCD on the lattice.To implement the exact chiral symmetry on the lattice, we use the optimal domain-wallfermion [8], of which the lattice fermion operator can be written as[ D ( m q )] xx (cid:48) ; ss (cid:48) ( m q ) = ( ω s D w + 1) xx (cid:48) δ ss (cid:48) + ( ω s D w − xx (cid:48) L ss (cid:48) , where { ω s , s = 1 , · · · , N s } are the exact solutions such that the effective 4-dimensional latticeDirac operator possesses the optimal chiral symmetry for any finite N s . The indices x and3 (cid:48) denote the lattice sites on the 4-dimensional lattice, and s and s (cid:48) the indices in the fifthdimension, while the Dirac and color indices have been suppressed. Here D w is the standardWilson Dirac operator plus a negative parameter − m (0 < m <
2) ( m is usually calledthe domain-wall height),( D w ) xx (cid:48) = (4 − m ) − (cid:88) ˆ µ =1 (cid:2) (1 − γ µ ) U µ ( x ) δ x +ˆ µ,x (cid:48) + (1 + γ µ ) U † µ ( x (cid:48) ) δ x − ˆ µ,x (cid:48) (cid:3) , where U µ ( x ) denotes the link variable pointing from x to x + ˆ µ . The operator L is independentof the gauge field, and it can be written as L = P + L + + P − L − , P ± = (1 ± γ ) / , (1)and ( L + ) ss (cid:48) = ( L − ) s (cid:48) s = − ( m q /m P V ) δ N s ,s (cid:48) , s = 1 ,δ s − ,s (cid:48) , < s ≤ N s , (2)where m q is the bare quark mass, and m P V = 2 m is the Pauli-Villars mass for the optimalDWF. Note that the matrices L ± satisfy L T ± = L ∓ , and R L ± R = L ∓ , where R is thereflection operator in the fifth dimension, with elements ( R ) ss (cid:48) = δ s (cid:48) ,N s +1 − s . Thus R L ± isreal and symmetric.Then the pseudofermion action for the optimal DWF can be written as S = φ † D ( m P V ) D ( m q ) φ, m P V = 2 m , where φ and φ † are complex scalar fields carrying the same quantum numbers (color, spin)of the fermion fields. Integrating the pseudofermion fields in the fermionic partition functiongives the fermion determinant of the effective 4-dimensional lattice Dirac operator D N s ( m q ),i.e., (cid:90) [ dφ † ][ dφ ] exp (cid:26) − φ † D ( m P V ) D ( m q ) φ (cid:27) = det D ( m q ) D ( m P V ) = det D N s ( m q ) , where D N s ( m q ) = m q + 12 ( m P V − m q )[1 + γ S N s ( H w )] ,S N s ( H w ) = 1 − (cid:81) N s s =1 T s (cid:81) N s s =1 T s , T s = 1 − ω s H w ω s H w .
4n the limit N s → ∞ , S N s ( H w ) → H w / (cid:112) H w , and D N s ( m q ) goes to D ( m q ) = m q + 12 ( m P V − m q ) (cid:34) γ H w (cid:112) H w (cid:35) . In the massless limit m q = 0, D (0) is equal to the overlap-Dirac operator [6], and it satisfiesthe Ginsparg-Wilson relation [9] D (0) γ + γ D (0) = 2 m P V D (0) γ D (0) ⇐⇒ D − γ + γ D − = 2 m P V γ , (3)where the chiral symmetry is broken by a contact term, i.e., the exact chiral symmetry atfinite lattice spacing. Note that (3) does not guarantee that any Ginsparg-Wilson Diracoperator D must possess exact zero modes in topologically non-trivial gauge background,not to mention to satisfy the Atiyah-Singer index theorem, Q t = n + − n − , where Q t is thetopological charge of the gauge background, and n ± is the number of exact zero modes of D with ± chirality. For example, the lattice Dirac operator constructed in Ref. [10] satisfiesthe Ginsparg-Wilson relation and possesses the correct axial anomaly in the continuum limit[11], but its index is always zero in any gauge background. So far, the overlap Dirac operatoris the only lattice Dirac operator to possesss topologically exact zero modes satisfying theAtiyah-Singer index theorem on a finite lattice.However, to perform HMC simulation of lattice QCD with the overlap Dirac operator isprohibitively expensive even for a small lattice (e.g., 16 × n ± at each step of the molecular dynamics [12].Moreover, the discontinuity of the fermion determinant at the topological boundary highlysuppresses the crossing rate between different topological sectors, thus renders HMC failingto sample all topological sectors ergodically. These difficulties can be circumvented by usingDWF with finite N s . Firstly, any positive lattice Dirac operator satisfying γ -hermiticity( γ Dγ = D † ) possesses a positive-definite pseudofermion action, without explicit depen-dence on n ± . Secondly, the step function of the fermion determinant at the topologicalboundary can be smoothed out by using DWF with finite N s (e.g., N s = 16), then theHMC on the 5-dimensional lattice can sample all topological sectors ergodically and alsokeep the chiral symmetry to a high precision with the optimal DWF [8, 13]. This has beendemonstrated for N f = 2 [14], N f = 1 + 1 [4], N f = 2 + 1 + 1 [15], and also N f = 2 + 1 + 1lattice QCD at the physical point [16]. 5 . Domain-wall fermion for heavy and light quarks In this subsection, we discuss which variant of DWF is more capable in capturing thequantum fluctuations of both heavy and light quarks in lattice QCD.Unlike other lattice fermions, DWF has the mass cutoff, i.e., the Pauli-Villars mass m P V ,and any quark mass has to satisfy the constraint m q (cid:28) m P V . Otherwise, if m q ∼ m P V , thendet( m q ) / det( m P V ) ∼
1, the internal quark loops are highly suppressed, and the quantumfluctuations of the quark field become mostly quenched. In general, the Pauli-Villars mass isequal to m P V a = 2 m (1 − dm ), where d is a parameter depending on the variant of DWF.For the Shamir[17]/M¨obius[18] DWF, d = 1 / m P V a = m (2 − m ) <
1, since m has tobe greater than 1 ( ∼ . − .
8) in order for its effective 4-dimensional Dirac operator to beable to detect the topology of a gauge configuration with nonzero topological charge. Thisimposes an upper-bound on the mass of Shamir/M¨obius heavy quark on the lattice, whichis more severe than the common constraint m q a < d = 0 and m P V a = 2 m (cid:29)
1, thusprovides the highest ceiling for accommodating the heavy quarks on the lattice, as well asthe minimal lattice artifacts due to the mass cutoff. This can be seen by comparing theeigenvalues of their effective 4D Dirac operators in the limit N s → ∞ , which is exactly equalto the overlap Dirac operator with the kernel H = cH w (1 + dγ H w ) − in the sign function, D ( m q ) = m q + 12 ( m P V − m q ) (cid:18) γ H √ H (cid:19) , m P V = 2 m (1 − dm ) , (4)where c = d = 1 / c = 1 and d = 0 for theBorici/Optimal DWF. The eigenvalues of (4) are lying on a circle in the complex planewith radius ( m P V − m q ) /
2, and center at m q + ( m P V − m q ) / m = 1 .
3, then m P V a = 2 m = 2 . m P V a = m (2 − m ) = 0 .
91 for the Shamir/M¨obius DWF. In Fig. 2, the eigenvaluesof (4) are plotted for m q = 0 (left panel) and m q a = 0 . m P V − m q ) / m q a = 0 . m P V − m q ) a = 0 .
11, and it shrinks tozero in the limit m q a → m P V a = 0 .
91. On the other hand, the Borici/Optimal DWF has m P V a = 2 m = 2 .
6, and ( m P V − m q ) a > m q a <
1, thus the eigenvalues of D ( m q )are not restricted to a very small circle even for the heavy quark. Moreover, in the chiral limit6 .0 0.5 1.0 1.5 2.0 2.5-1.5-1.0-0.50.00.51.01.5 Borici / Optimal DWFShamir / Mobius DWF
Re( ) I m ( ) m q = 0.0m = 1.3 Borici / Optimal DWFShamir / Mobius DWF
Re( ) I m ( ) m q a = 0.8m = 1.3 FIG. 2: Comparing the eigenvalue spectra of the effective 4D Dirac operator of the Borici/OptimalDWF and the Shamir/M¨obius DWF for m q = 0 (left panel) and m q a = 0 . (left panel), the radius of the eigenvalue circle for the Borici/Optimal DWF is more than2 times of that of the Shamir/M¨obius DWF. This implies that the Borici/Optimal DWFis more capable than the Shamir/M¨obius DWF in capturing the short-distance quantumfluctuations of the QCD vacuum, for both light and heavy quarks. C. Zolotarev optimal rational approximation and optimal domain-wall fermion
For any numerical simulation of lattice QCD with DWF, an important question is whatis the optimal chiral symmetry for any finite N s in the fifth dimension, in the sense how itseffective 4D lattice Dirac operator can be exactly equal to the Zolotarev optimal rationalapproximation of the overlap Dirac operator. The exact solution to this problem is given inRef. [8], with the optimal { ω s } ω s = 1 λ min (cid:112) − κ (cid:48) sn ( v s ; κ (cid:48) ) , s = 1 , · · · , N s , (5)where sn( v s ; κ (cid:48) ) is the Jacobian elliptic function with argument v s (see Eq. (13) in Ref.[8]) and modulus κ (cid:48) = (cid:112) − λ min /λ max . Then S N s ( H w ) is exactly equal to the Zolotarevoptimal rational approximation of H w / (cid:112) H w , i.e., the approximate sign function S N s ( H w )7atisfying the bound | − S N s ( λ ) | ≤ d Z for λ ∈ [ λ min , λ max ], where d Z is the maximumdeviation | − √ xR Z ( x ) | max of the Zolotarev optimal rational polynomial R Z ( x ) of 1 / √ x for x ∈ [1 , λ max /λ min ], with degree ( n − , n ) for N s = 2 n .Nevertheless, the optimal weights { ω s } in (5) do not satisfy the R symmetry ( ω s = ω N s − s +1 ) which is required for the exact one-flavor pseudofermion action for DWF [4]. Theoptimal { ω s } satisfying R symmetry is obtained in Ref. [13]. For N s = 2 n , the optimal { ω s } satisfying R symmetry are written as ω s = ω N s +1 − s = 1 λ min (cid:115) − κ (cid:48) sn (cid:18) (2 s − K (cid:48) N s ; κ (cid:48) (cid:19) , s = 1 , · · · , N s / , (6)where sn( u ; κ (cid:48) ) is the Jacobian elliptic function with modulus κ (cid:48) = (cid:112) − λ min /λ max , and K (cid:48) is the complete elliptic function of the first kind with modulus κ (cid:48) . Then the approximatesign function S N s ( H w ) satisfies the bound 0 ≤ − S N s ( λ ) ≤ d Z for λ ∈ [ λ min , λ max ], where d Z is defined above. Note that δ ( λ ) = 1 − S ( λ ) does not satisfy the criterion that the maximaand minima of δ ( λ ) all have the same magnitude but with the opposite sign ( δ min = − δ max ).However, the most salient featues of the optimal rational approximation of degree ( m, n ) arepreserved, namely, the number of alternate maxima and minima is ( m + n + 2), with ( n + 1)maxima and ( m + 1) minima, and all maxima (minma) are equal to 2 d Z (0). This can beregarded as the generalized optimal rational approximation (with a constant shift).In this study, the parameters for the pseudofermion action are: m = 1 . N s = 2 n = 16, λ max /λ min = 6 . / .
05, and the optimal weights { ω s , s = 1 , · · · , N s } for the 2-flavor partsare obtained with (5), while for the one-flavor parts with (6). In Fig. 3, the deviation ofthe sign function, δ ( λ ) = 1 − S ( λ ), is plotted versus λ , for (a) without the R symmetry,and (b) with the R symmetry. Here δ ( λ ) has 2 n + 1 = 17 alternate maxima and minimain the interval [ λ min , λ max ] = [0 . , . − d Z ≤ − S ( λ ) ≤ d Z , while for (b), 0 ≤ − S ( λ ) ≤ d Z , where d Z is the maximum deviation | − √ xR (7 , Z | max of the Zolotarev optimal rational polynomial. III. GENERATION OF THE GAUGE ENSEMBLE
In this section, we give the details of the actions, the algorithms, and the parametersto perform the HMC simulations in this study. Moreover, for the initial 257 trajectoriesgenerated by a single node (with 2 Nvidia GTX-TITAN-X GPU cards), the topological8 ( ) = S ( ) -15x10 -6 -10x10 -6 -5x10 -6 -6 -6 -6 ( ) = S ( ) -6 -6 -6 -6 -6 (a) (b)FIG. 3: The deviation δ ( λ ) = 1 − S ( λ ) of the optimal DWF with N s = 2 n = 16 and λ max /λ min =6 . / .
05, for (a) without R symmetry, and (b) with R symmetry. charge fluctuation is measured, and the HMC characteristics are presented. Details of thelattice setup are given as follows. A. The actions
In the following, we present the details of the fermion actions and the gauge action inour HMC simulations.As noted in Ref. [15], for domain-wall fermions (DWF), to simulate N f = 2 + 1 + 1amounts to simulate N f = 2 + 2 + 1. Similarly, to simulate N f = 2 + 1 + 1 + 1 amounts tosimulate N f = 2 + 2 + 1 + 1, i.e., (cid:18) det D ( m u/d )det D ( m P V ) (cid:19) det D ( m s )det D ( m P V ) det D ( m c )det D ( m P V ) det D ( m b )det D ( m P V )= (cid:18) det D ( m u/d )det D ( m P V ) (cid:19) (cid:18) det D ( m c )det D ( m P V ) (cid:19) det D ( m s )det D ( m c ) det D ( m b )det D ( m P V ) , (7)where only one of the 6 possible possibilities for N f = 2+2+1+1 is written. Note that on theRHS of Eq. (7), the 2-flavor simulation with (det D ( m c ) / det D ( m P V )) is more efficient thanits counterpart of one-flavor with (det D ( m c ) / det D ( m P V )) on the LHS. Moreover, the one-flavor simulation with det D ( m s ) / det D ( m c ) on the RHS is more efficient than the original9ne with det D ( m s ) / det D ( m P V ) on the LHS. Thus, we perform the HMC simulation withthe expression on the RHS of Eq. (7).For the two-flavor parts, (cid:0) det D ( m u/d ) / det D ( m P V ) (cid:1) and (det D ( m c ) / det D ( m P V )) , weuse the N f = 2 pseudofermion action which has been using since 2011 [14], and it can bewritten as S ( m q , m P V ) = φ † C † ( m P V ) { C ( m q ) C † ( m q ) } − C ( m P V ) φ, m P V = 2 m , (8)where C ( m q ) = 1 − M ( m q ) D OE w M ( m q ) D EO w ,M ( m q ) = { − m + ω − / [1 − L ( m q )][(1 + L ( m q )] − ω − / } − , and L ( m q ) is defined in (1) and (2). Here ω ≡ diag { ω , ω , · · · , ω N s } is a diagonal matrixin the fifth dimension, and D EO / OE w denotes the part of D w with gauge links pointing fromeven/odd sites to odd/even sites after even-odd preconditioning on the 4-dimensional lattice.For the two-flavor part of u and d quarks, we turn on the mass-preconditioning [20]by introducing an auxillary heavy fermion field with mass m H a = 0 .
1. Then the N f = 2pseudofermion action (8) is replaced with S ( m q , m H ) + S ( m H , m P V )= φ † C ( m H ) † { C ( m q ) C ( m q ) † } − C ( m H ) φ + φ † H C † ( m P V ) { C ( m H ) C ( m H ) † } − C ( m P V ) φ H , which gives the partition function (fermion determinant) exactly the same as that of (8).For the one-flavor parts, det D ( m s ) / det D ( m c ) and det D ( m b ) / det D ( m P V ), we use theexact one-flavor pseudofermion action (EOFA) for DWF [4]. For the optimal DWF, it canbe written as ( m < m )det D ( m )det D ( m ) = det D T ( m )det D T ( m ) = (cid:90) dφ †± dφ ± exp (cid:16) − φ † + G + ( m , m ) φ + − φ †− G − ( m , m ) φ − (cid:17) , (9)where φ ± and φ †± are pseudofermion fields (each of two spinor components) on the 4-dimensional lattice, and G − ( m , m ) = P − (cid:20) I − k ( m , m ) ω − / v T − H T ( m ) v − ω − / (cid:21) P − , (10) G + ( m , m ) = P + (cid:20) I + k ( m , m ) ω − / v T + H T ( m ) − ∆ + ( m , m ) P + v + ω − / (cid:21) P + . (11)10ere D T ( m i ) = D w + M ( m i ) , i = 1 , M ( m i ) = ω − / [1 − L ( m i )][1 + L ( m i )] − ω − / = P + M + ( m i ) + P − M − ( m i ) ,H T ( m i ) = R γ D T ( m i ) , ∆( m , m ) = R [ M ( m ) − M ( m )] = P + ∆ + ( m , m ) + P − ∆ − ( m , m ) , ∆ ± ( m , m ) = k ( m , m ) ω − / v ± v T ± ω − / ,k ( m , m ) = m − m m + m ,v T + = ( − , , · · · , ( − N s ) , v − = − v + . For the gluon fields, we use the Wilson plaquette gauge action [22] at β = 6 /g = 6 . S g ( U ) = 6 g (cid:88) plaq. (cid:26) −
13 ReTr( U p ) (cid:27) , where g is the bare coupling.The bare mass of u / d quarks is set to m u/d = 0 .
01 such that M π L >
4, while the baremasses of ( b , c , s ) are tuned to ( m b , m c , m s ) = (0 . , . , .
15) such that they give themasses of the vector mesons Υ(9460),
J/ψ (3097) , and φ (1020) respectively. The algorithmfor simulating 2-flavor action for optimal domain-wall quarks has been outlined in Ref. [14],while that for simulating the exact one-flavor pseudofermion action (EOFA) of domain-wallfermion has been presented in Refs. [4, 21]. In the molecular dynamics, we use the Omelyanintegrator [23], the multiple-time scale method [24], and the mass-preconditioning [20]. B. HMC simulations
Following the common strategy to reduce the thermalization time for a large lattice suchas 40 ×
64, we first perform the thermalization on a smaller lattice 20 ×
32 with thesame set of parameters ( β, m u/d , m s , m c , m b ). Then the thermalized gauge configurationon the 20 ×
32 lattice is used to construct the inital gauge configuration on the 40 × × rajectory number
50 100 150 200 250 - F l a v o r F o r c e s - F l a v o r F o r c e s G auge F o r c e (m H / m PV ) (m c /m PV ) (m u /m H ) (m s /m c ) (m s /m c ) (m b /m PV ) (m b /m PV ) FIG. 4: The maximum forces of the gauge field, the 2-flavor pseudofermion fields, and the one-flavor pseudofermion fields versus the HMC trajectory in the HMC simulations of the lattice QCDwith N f = 2 + 1 + 1 + 1 optimal DWF. trajectories, resulting 14 “seed” configurations. Then we use these seed configurations as theinitial configurations for 14 independent simulations on 14 nodes, each of two Nvidia GTX-TITAN-X GPU cards. Each node generates ∼
40 trajectories independently, and all 14 nodesaccumulate a total of 535 trajectories. We sample one configuration every 5 trajectories ineach stream, and obtain a total of 103 configurations for physical measurements.In the following, we summarize the HMC characteristics of the first 257 trajectories. InFig. 4, we plot the maximum force (averaged over all links) among all momentum updatesin each trajectory, for the gauge force, the 2-flavor pseudofermion forces, and the one-flavor pseudofermion forces respectively, where φ ( m /m ) denotes the two-flavor fermionforce due to the pseudofermion action S ( m , m ), and φ ± ( m /m ) denotes the one-flavorpseudofermion force due to the exact one-flavor action with ± chirality, S ± ( m , m ) =12 arjectory number
50 100 150 200 250 H -4-2024 FIG. 5: The change of the Hamiltonian ∆ H versus the trajectory in the HMC simulations of latticeQCD with N f = 2 + 1 + 1 + 1 optimal DWF. The line connecting the data points is only for guidingthe eyes. φ †± G ± ( m , m ) φ ± . From the sizes of various forces in Fig. 4, the multiple time scales canbe designed in the momentum update with the gauge force and the pseudofermion forces.With the length of the HMC trajectory equal to one, we use 4 different time scales for themomentum updates with (1) the gauge force; (2) the two-flavor fermion forces associatedwith φ ( m c /m P V ) and φ ( m H /m P V ); (3) the two-flavor force associated with φ ( m u /m H ) andthe one-flavor fermion force associated with φ + ( m b /m P V ); (4) the one-flavor fermion forcesassociated with φ − ( m b /m P V ), φ + ( m s /m c ), and φ − ( m s /m c ), which correspond to the stepsizes 1 / ( k k k k ), 1 / ( k k k ), 1 / ( k k ), and 1 /k respectively. In our simulation, we set( k , k , k , k ) = (10 , , , H versus the HMC trajectory is plotted for the first257 trajectories, with (cid:104) ∆ H (cid:105) = 0 . . (cid:104) ∆ H (cid:105) = 0 . P acc = erfc (cid:16)(cid:112) (cid:104) ∆ H(cid:105) / (cid:17) [25], which gives 0.664(24), in good agreement with the measured acceptance rate 0.673(29).Moreover, we measure the expectation value of exp( − ∆ H ), to check whether it is consistentwith the theoretical formula (cid:104) exp( − ∆ H ) (cid:105) = 1 which follows from the area-preserving prop-erty of the HMC simulation [26]. The measured value of (cid:104) exp( − ∆ H ) (cid:105) is 1.026(66), in good13greement with the theoretical expectation value. The summary of the HMC characteristicsfor the inital 257 trajectories is given in Table I. TABLE I: Summary of the HMC characteristics for the first 257 trajectories in the simulation of N f = 2 + 1 + 1 + 1 lattice QCD with the optimal DWF. N traj Time(s)/traj Acceptance (cid:104) ∆ H (cid:105) P acc = erfc( (cid:112) (cid:104) ∆ H (cid:105) / (cid:104) exp( − ∆ H ) (cid:105) (cid:104) plaquette (cid:105)
257 76349(146) 0.673(29) 0.376(57) 0.664(24) 1.026(66) 0.63185(1)
C. Topological charge fluctuations
In this subsection, we examine the evolution of the topological charge Q t in the first 257trajectories, and obtain the histogram of its distribution.In lattice QCD with exact chiral symmetry, the topological charge Q t can be measured bythe index of the massless overlap-Dirac operator, since its index satisfies the Atiyah-Singerindex theorem, index( D ov ) = n + − n − = Q t . However, to project the zero modes of themassless overlap-Dirac operator for the 40 ×
64 lattice is prohibitively expensive. On theother hand, the clover topological charge Q clover = (cid:80) x (cid:15) µνλσ tr[ F µν ( x ) F λσ ( x )] / (32 π ) is notreliable [where the matrix-valued field tensor F µν ( x ) is obtained from the four plaquettessurrounding x on the (ˆ µ, ˆ ν ) plane], unless the gauge configuration is sufficiently smooth.Nevertheless, the smoothness of a gauge configuration can be attained by the Wilson flow[27, 28], which is a continuous-smearing process to average gauge field over a sphericalregion of root-mean-square radius R rms = √ t , where t is the flow-time. In this study, theflow equation is numerically integrated from t = 0 with ∆ t/a = 0 .
01, and measure the Q clover at t/a = 0 . R rms = √ t ∼ . a . Then each gauge configuration becomesvery smooth, with Q clover close to an integer, and the average plaquette greater than 0.997.Denoting the nearest integer of Q clover by Q t ≡ round( Q clover ), Q t is plotted versus thetrajectory number in the left-panel of Fig. 6, while the right-panel displays the histogramof the probability distribution of Q t of the first 257 HMC trajectories. Evidently, the HMCsimulation samples all topological sectors ergodically. Note that the lattice volume is rather14mall, ∼ (1.2 fm) × (1.9 fm), thus the chance of sampling Q t > trajectory
50 100 150 200 250 Q t -4-2024 Q t -4 -3 -2 -1 0 1 2 3 4 P ( Q t ) FIG. 6: (left panel) The evolution of Q t versus the HMC trajectory. The line connecting the datapoints is only for guiding the eyes. (right panel) The histogram of the probability distribution of Q t for the first 257 HMC trajectories. D. Lattice scale
First, we recap the generation of the gauge ensemble. From the initial 257 trajectoriesgenerated by a single node, we discard the first 187 trajectories for thermalization, andsample one configuration every 5 trajectories, resulting 14 “seed” configurations. Then weuse these seed configurations as the initial configurations for 14 independent simulationson 14 nodes, each of two Nvidia GTX-TITAN-X GPU cards. Each node generates ∼ { t (cid:104) E ( t ) (cid:105)} (cid:12)(cid:12) t = t = 0 . , and obtain √ t /a = 4 . √ t = 0 . a − = 6 . ± .
037 GeV. The lattice spacing is a = 0.0303(2) fm, giving thespatial volume ∼ (1.213 fm) , which is too small for studying physical observables involvingthe light quarks. E. Quark propagator
We compute the valence quark propagator of the effective 4D Dirac operator with thepoint source at the origin, and with the mass and other parameters exactly the same asthose of the sea quarks. The boundary conditions are periodic in space and antiperiodic intime. First, we solve the following linear system with mixed-precision conjugate gradientalgorithm, for the even-odd preconditioned D [30] D ( m q ) | Y (cid:105) = D ( m P V ) B − | source vector (cid:105) , (12)where B − x,s ; x (cid:48) ,s (cid:48) = δ x,x (cid:48) ( P − δ s,s (cid:48) + P + δ s +1 ,s (cid:48) ) with periodic boundary conditions in the fifthdimension. Then the solution of (12) gives the valence quark propagator( D c + m q ) − x,x (cid:48) = ( m P V − m q ) − [( BY ) x, x (cid:48) , − δ x,x (cid:48) ] , m P V = 2 m . (13)Each column of the quark propagator is computed by a single node with 2 Nvidia GTX-TITAN-X GPU cards, which attains more than 1000 Gflops/sec (sustained). F. Residual masses
To measure the chiral symmetry breaking due to finite N s , we compute the residual massaccording to [31], m res = (cid:42) tr( D c + m q ) − , tr[ γ ( D c + m q ) γ ( D c + m q )] − , (cid:43) − m q , (14)where ( D c + m q ) − denotes the valence quark propagator with m q equal to the sea-quarkmass, tr denotes the trace running over the color and Dirac indices, and the brackets (cid:104)· · · (cid:105) denote the averaging over the gauge ensemble. In the limit N s → ∞ , D c is exactly chiralsymmetric and the first term on the RHS of (14) is exactly equal to m q , thus the residual mass m res is exactly zero, and the quark mass m q is well-defined for each gauge configuration.On the other hand, for any finite N s with nonzero residual mass, the quark mass is not16ell-defined for each gauge configuration, but its impact on any physical observable can beroughly estimated by the difference due to changing the valence quark mass from m q to m q + m res . TABLE II: The residual masses of u / d , s , c , and b quarks.quark m q a m res a m res [MeV] u / d . × − s . × − c . × − b . × − For the 103 gauge configurations generated by HMC simulation of lattice QCD with N f = 2 + 1 + 1 + 1 optimal domain-wall quarks, the residual masses of u / d , s , c , and b quarks are listed in Table II. We see that the residual mass of any quark flavor is less than0 .
007 MeV, which should be negligible in comparison with other systematic uncertainties.
IV. MASS SPECTRA OF BEAUTY MESONS
In the following, we determine the masses of the low-lying mesons with valence quarkcontents ¯bb , ¯bc , ¯bs , and ¯cc . We construct the quark-antiquark meson interpolators andmeasure their time-correlation functions using the point-to-point quark propagators com-puted with the same parameters ( N s = 16, m = 1 . λ max /λ min = 6 . / .
05) of the seaquarks, for the quark masses ( m u/d a = 0 . m s a = 0 . m c a = 0 . m b a = 0 . m b , m c and m s are fixed by the masses of the vector mesons Υ(9460), J/ψ (3097), and φ (1020) respectively. Then we extract the mass of the lowest-lying meson state from thetime-correlation function.The time-correlation function of the beauty meson interpolator ¯b Γ q (where q = { b , c , s } )is measured according to the formula C Γ ( t ) = (cid:42)(cid:88) (cid:126)x tr { Γ( D c + m b ) − x, Γ( D c + m q ) − ,x } (cid:43) , (15)where Γ = { , γ , γ i , γ γ i , (cid:15) ijk γ j γ k } , corresponding to scalar ( S ), pseudoscalar ( P ), vector( V ), axial-vector ( A ), and pseudovector ( T ) respectively, and the valence quark propagator17 D c + m q ) − is computed according to the formula (13). Note that ¯q γ γ i q transforms like J P C = 1 ++ , while ¯q (cid:15) ijk γ j γ k q like J P C = 1 + − .For the vector meson, we average over i = 1 , , C V ( t ) = (cid:42) (cid:88) i =1 (cid:88) (cid:126)x tr { γ i ( D c + m b ) − x, γ i ( D c + m q ) − ,x } (cid:43) . Similarly, we perform the same averaging for the axial-vector and pseudovector mesons.Moreover, to enhance statistics, we average the forward and the backward time-correlationfunction. ¯ C ( t ) = 12 [ C ( t ) + C ( T − t )] . The time-correlation function (TCF) and the effective mass of the meson interpolators ¯b Γ b , ¯c Γ c , ¯b Γ c , and ¯b Γ s are plotted in the Appendices A-D respectively. A. Bottomonium and Charmonium
First of all, we check to what extent we can reproduce the bottomonium masses whichhave been measured precisely in high energy experiments.Our results of the mass spectrum of the low-lying states of bottomonium are summarizedin Table III. The time-correlation function and the effective mass of ¯b Γ b are plotted inAppendix A.The first column in Table III is the Dirac matrix used for computing the time-correlationfunction (15). The second column is J P C of the state. The third column is the [ t , t ] usedfor fitting the data of C Γ ( t ) to the usual formula z M a [ e − Mat + e − Ma ( T − t ) ] (16)to extract the ground state meson mass M , where the excited states have been neglected.We use the correlated fit throughout this work. The fifth column is the mass M of themeson state, where the first error is statistical, and the second is systematic. Here thestatistical error is estimated using the jackknife method with the bin size of which thethe statistical error saturates, while the systematic error is estimated based on all fittingssatisfying χ / dof < . | t − t | ≥ t ≥
10 and t ≤
32. The last column is theexperimental state we have identified, and its PDG mass value [32].18he analysis and the descriptions in the above paragraph apply to all results obtained inthis work, as given in Table III-VI.
TABLE III: The masses of low-lying bottomonium states obtained in this work. The fifth columnis the mass of the meson state, where the first error is statistical, and the second is systematic.The last column is the experimental state we have identified, and its PDG mass value [32]. For adetailed description of each column, see the paragraph with Eq. (16).Γ J P C [ t , t ] χ /dof Mass(MeV) PDG1I 0 ++ [19,29] 1.10 9859(14)(11) χ b (9859) γ − + [15,31] 1.04 9403(4)(5) η b (9399) γ i −− [21,31] 0.51 9468(7)(6) Υ(9460) γ γ i ++ [19,26] 1.15 9884(27)(35) χ b (9893) (cid:15) ijk γ j γ k + − [19,25] 0.97 9910(20)(25) h b (9899) Evidently, the masses of bottomonium in Table III are in good agreement with the PDGmass values, even though the axial-vector (1 −− ) and pseudovector (1 + − ) mesons have rela-tively larger errors than other meson states. TABLE IV: The masses of low-lying charmonium states obtained in this work. The fifth columnis the mass of the meson state, where the first error is statistical, and the second is systematic.The last column is the experimental state we have identified, and its PDG mass value [32]. For adetailed description of each column, see the paragraph with Eq. (16).Γ J P C [ t , t ] χ /dof Mass(MeV) PDG1I 0 ++ [14,25] 1.01 3403(16)(13) χ c (3415) γ − + [15,29] 1.17 2989(6)(4) η c (2984) γ i −− [15,28] 0.65 3112(7)(5) J/ψ (3097) γ γ i ++ [14,21] 1.13 3513(23)(10) χ c (3510) (cid:15) ijk γ j γ k + − [17,25] 0.39 3527(14)(19) h c (3524) Next, we turn to the charmonium states extracted from the gound states of ¯c Γ c . Our19esults of the masses of the low-lying states of charmonium are summarized in Table IV.The time-correlation function and the effective mass of ¯c Γ c are plotted in Appendix B.Evidently, the theoretical masses of charmonium in Table IV are in good agreement withthe PDG values. Note that the theoretical result of the hyperfine splitting (1 S − S ) is123(9)(6) MeV, in good agreement with the PDG value 118 MeV. B. B s and B c mesons Our results of the masses of the low-lying states of B s mesons are summarized in TableV. The time-correlation function and the effective mass of ¯b Γ s are plotted in AppendixD. Here we have identified the scalar ¯bs meson with the state B ∗ sJ (5850) observed in highenergy experiments, due to the proximity of their masses. This predicts that B ∗ sJ (5850)possesses J P = 0 + , which can be verified by high energy experiments in the future. Moreover,the pseudovector meson (the last entry in Table V) has not been observed in high energyexperiments, thus it serves as a prediction of N f = 2 + 1 + 1 + 1 lattice QCD. TABLE V: The masses of low-lying B s meson states obtained in this work. The fifth column is themass of the meson state, where the first error is statistical, and the second is systematic. The lastcolumn is the experimental state we have identified, and its PDG mass value [32]. For a detaileddescription of each column, see the paragraph with Eq. (16).Γ J P [ t , t ] χ /dof Mass(MeV) PDG1I 0 + [15,24] 0.37 5839(30)(18) B ∗ sJ (5850) γ − [23,29] 0.79 5406(16)(17) B s (5367) γ i − [18,29] 0.66 5430(17)(18) B ∗ s (5415) γ γ i + [16,22] 0.58 5839(23)(14) B s (5830) (cid:15) ijk γ j γ k + [16,23] 0.56 5909(26)(34) Finally, we turn to the heavy mesons with beauty and charm. In Table VI, we summarizeour results of the masses of B c mesons extracted from the ground states of ¯b Γ c . The time-correlation function and the effective mass of ¯b Γ c are plotted in the Appendix C. Exceptfor the pseudoscalar meson B c (6275), other four meson states have not been observed in20igh energy experiments. It is interesting to see to what extent the experimental results willagree with our theoretical predictions. TABLE VI: The masses of low-lying B c meson states obtained in this work. The fifth column is themass of the meson state, where the first error is statistical, and the second is systematic. The lastcolumn is the experimental state we have identified, and its PDG mass value [32]. For a detaileddescription of each column, see the paragraph with Eq. (16).Γ J P [ t , t ] χ /dof Mass(MeV) PDG1I 0 + [20,28] 1.17 6766(38)(16) γ − [15,31] 1.02 6285(6)(5) B c (6275) γ i − [16,31] 0.68 6375(6)(7) γ γ i + [21,32] 0.62 6787(34)(28) (cid:15) ijk γ j γ k + [19,26] 0.97 6798(33)(17) V. QUARK MASSES OF (b , c , s) The quark masses cannot be measured directly in high energy experiments since quarksare confined inside hadrons. Therefore, the quark masses can only be determined by com-paring theoretical calculations of physical observables with the experimental values. Forany field theoretic calculation, the quark masses depend on the regularization, as well asthe renormalization scheme and scale. For lattice QCD, the hadron masses can be com-puted nonperturbatively from the first principles, and from which the quark masses can bedetermined.We have used the mass of the vector meson Υ(9460) to fix the bare mass of b quark equalto m b = 0 . a − . To transcribe the bare mass to the corresponding value in the usualrenormalization scheme MS in high energy phenomenology, one needs to compute the latticerenormalization constant Z m = Z − s , where Z s is the renormalization constant for ¯ ψψ . Ingeneral, Z m should be determined nonperturbatively. However, in this study, the latticespacing is rather small ( a (cid:39) .
03 fm), thus it is justified to use the one-loop perturbation21ormula [33] Z s ( µ ) = 1 + g π (cid:2) ln( a µ ) + 0 . (cid:3) ( m = 1 . . (17)At β = 6 . a − = 6 . µ = 2 GeV, (17) gives Z s = 1 . m b to the MS mass at µ = 2 GeV m b (2 GeV) = m b Z m (2 GeV) = 5 . ± .
025 GeV , where the error bar combines (in quadrature) the statistical error and the systematic errorsof the lattice spacing and the b quark bare mass.To compare our result with the PDG value of m b ( m b ) at the scale µ = m b , we solve theequation m b = m b Z m ( µ = m b ) and obtain m b ( m b ) = 4 . ± .
04 GeV , (18)which is higher than the PDG value (4 . ± .
03) GeV for N f = 2 + 1 + 1 lattice QCD, butis closer to the value in the 1S scheme m b = 4 . m c =0 . a − is transcribed to m c (2 GeV) = 1 . ± .
03 GeV , where the error bar combines (in quadrature) the statistical and the systematic errors fromthe lattice spacing and the charm quark bare mass. To compare our result with the PDGvalue of m c ( m c ), we solve m c = m c Z m ( µ = m c ) and obtain m c ( m c ) = 1 . ± .
03 GeV , (19)which is slightly smaller than the PDG value (1 . ± . N f = 2 + 1 + 1 latticeQCD [32].Finally we turn to the strange quark mass. Using (17), the strange quark bare mass m s = 0 . a − is transcribed to m s (2 GeV) = 88 . ± . , (20)where the error bar combines (in quadrature) the statistical and the systematic ones fromthe lattice spacing and the s quark bare mass. Our result of the strange quark mass (20) isslightly smaller than the PDG value (92 . ± .
7) MeV for N f = 2 + 1 + 1 lattice QCD [32].22 I. CONCLUDING REMARK
This study demonstrates that the Dirac b quark can be simulated dynamically in latticeQCD, together with the ( c , s , d , u ) quarks. Even with unphysically heavy u and d quarks inthe sea, the low-lying mass spectra of mesons with valence quark contents ¯bb , ¯bc , ¯bs , and ¯cc are in good agreement with the experimental values. Also, we have several predictionswhich have not been observed in high energy experiments, i.e., predicting the mass and the J P of four B c meson states (see Table VI), the J P of B ∗ sJ (5850) to be 0 + , and the mass andthe J P of the pseudovector B s meson state (see Table V). Moreover, we have determinedthe masses of ( b , c , s ) quarks, as given in (18), (19), and (20) respectively.These results imply that it is feasible to simulate lattice QCD with physical ( u , d , s , c , b )domain-wall quarks on a large ( ∼ ) lattice, with the Exaflops supercomputers which willbe available ∼ u , d , s , c , b ) quark contents can becomputed from the first principles of QCD. This will provide a viable way to systematicallyreduce the uncertainties in the theoretical predictions of the Standard Model (SM), whichare largely stemming from the sector of the strong interaction[34]. This is crucial for un-veiling any new physics beyond the standard model (SM), by identifying any discrepanciesbetween the high energy experimental results and the theoretical values derived from thefirst principles of the SM with all quarks (heavy and light) as Dirac fermions, without usingnon-relativistic approximation or heavy quark effective field theory for b and c quarks. Acknowledgements
The author is grateful to Academia Sinica Grid Computing Center (ASGC) and NationalCenter for High Performance Computing (NCHC) for the computer time and facilities. Thiswork is supported by the Ministry of Science and Technology (Grant Nos. 108-2112-M-003-005, and 107-2119-M-003-008). [1] T. W. Chiu et al. [TWQCD Collaboration], Phys. Lett. B , 171 (2007) [arXiv:0705.2797[hep-lat]];
2] T. W. Chiu et al. [TWQCD Collaboration], PoS LAT , 180 (2007) [arXiv:0704.3495[hep-lat]].[3] S. Duane, A. D. Kennedy, B. J. Pendleton and D. Roweth, Phys. Lett. B , 216 (1987).[4] Y. C. Chen, T. W. Chiu [TWQCD Collaboration], Phys. Lett. B , 55 (2014)[arXiv:1403.1683 [hep-lat]].[5] D. B. Kaplan, Phys. Lett. B , 342 (1992) [hep-lat/9206013].[6] H. Neuberger, Phys. Lett. B , 141 (1998) [hep-lat/9707022].[7] R. Narayanan and H. Neuberger, Nucl. Phys. B , 305 (1995) [hep-th/9411108].[8] T. W. Chiu, Phys. Rev. Lett. , 071601 (2003) [hep-lat/0209153][9] P. H. Ginsparg and K. G. Wilson, Phys. Rev. D , 2649 (1982).[10] T. W. Chiu, Phys. Lett. B , 429 (2001) [hep-lat/0106012].[11] T. W. Chiu and T. H. Hsieh, Phys. Rev. D , 054508 (2002) [hep-lat/0109016].[12] Z. Fodor, S. D. Katz and K. K. Szabo, JHEP , 003 (2004)[13] T. W. Chiu, Phys. Lett. B , 95 (2015) [arXiv:1503.01750 [hep-lat]].[14] T. W. Chiu, T. H. Hsieh, Y. Y. Mao [TWQCD Collaboration], Phys. Lett. B , 420 (2012)[arXiv:1109.3675 [hep-lat]].[15] Y. C. Chen and T. W. Chiu [TWQCD Collaboration], Phys. Lett. B , 193 (2017)[arXiv:1701.02581 [hep-lat]].[16] T. W. Chiu [TWQCD Collaboration], arXiv:2002.06126 [hep-lat].[17] Y. Shamir, Nucl. Phys. B , 90 (1993) [hep-lat/9303005].[18] R. C. Brower, H. Neff and K. Orginos, Nucl. Phys. Proc. Suppl. , 686 (2005) [hep-lat/0409118].[19] A. Borici, Nucl. Phys. Proc. Suppl. (2000) 771 [hep-lat/9909057].[20] M. Hasenbusch, Phys. Lett. B , 177 (2001) [hep-lat/0107019].[21] Y. C. Chen, T. W. Chiu [TWQCD Collaboration], PoS IWCSE , 059 (2014)[arXiv:1412.0819 [hep-lat]].[22] K. G. Wilson, Phys. Rev. D , 2445 (1974).[23] I. P. Omelyan, I. M. Mryglod, and R. Folk, Phys. Rev. Lett. , 898 (2001).[24] J. C. Sexton and D. H. Weingarten, Nucl. Phys. B , 665 (1992).[25] S. Gupta, A. Irback, F. Karsch and B. Petersson, Phys. Lett. B , 437 (1990).[26] M. Creutz, Phys. Rev. D , 1228 (1988).
27] R. Narayanan and H. Neuberger, JHEP , 064 (2006) [hep-th/0601210].[28] M. Luscher, JHEP , 071 (2010) Erratum: [JHEP , 092 (2014)] [arXiv:1006.4518[hep-lat]].[29] A. Bazavov et al. [MILC Collaboration], Phys. Rev. D , no. 9, 094510 (2016)[arXiv:1503.02769 [hep-lat]].[30] T. W. Chiu et al. [TWQCD Collaboration], PoS LATTICE , 030 (2010) [arXiv:1101.0423[hep-lat]].[31] Y. C. Chen, T. W. Chiu [TWQCD Collaboration], Phys. Rev. D , 094508 (2012)[arXiv:1205.6151 [hep-lat]].[32] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D , no. 3, 030001 (2018) and 2019update.[33] C. Alexandrou, E. Follana, H. Panagopoulos and E. Vicari, Nucl. Phys. B , 394 (2000).[34] The t quark can be neglected in the strong interaction since it is very short-lived and it decaysto W -boson and b / s / d quarks before it can interact with other quarks through the gluons. ppendix A: C ( t ) and the effective mass of ¯b Γ b t C ( t ) -5 -10 -15 -20 -25 (pseudoscalar) _b b t M e ff a (pseudoscalar) _b b FIG. 7: The time-correlation function and the effective mass of the meson interpolator ¯b γ b . t C ( t ) -5 -10 -15 -20 -25 i (vector) _b b t M e ff a i (vector) _b b FIG. 8: The time-correlation function and the effective mass of the meson interpolator ¯b γ i b . C ( t ) -5 -10 -15 -20 -25 (scalar) _b b t M e ff a (scalar)_b b FIG. 9: The time-correlation function and the effective mass of the meson interpolator ¯bb . t C ( t ) -5 -10 -15 -20 -25 i (axial-vector) _b b t M e ff a i (axial-vector)_b b FIG. 10: The time-correlation function and the effective mass of the meson interpolator ¯b γ γ i b . C ( t ) -5 -10 -15 -20 -25 ijk j k (tensor) _b b t M e ff a ijk j k (tensor)_b b FIG. 11: The time-correlation function and the effective mass of the meson interpolator ¯b (cid:15) ijk γ j γ k b . ppendix B: C ( t ) and the effective mass of ¯c Γ c t C ( t ) -2 -4 -6 -8 -10 -12 (pseudoscalar) _c c t M e ff a (pseudoscalar)_c c FIG. 12: The time-correlation function and the effective mass of the meson interpolator ¯c γ c . t C ( t ) -2 -4 -6 -8 -10 -12 i (vector) _c c t M e ff a i (vector)_c c FIG. 13: The time-correlation function and the effective mass of the meson interpolator ¯c γ i c . C ( t ) -2 -4 -6 -8 -10 -12 (scalar) _c c t M e ff a (scalar)_c c FIG. 14: The time-correlation function and the effective mass of the meson interpolator ¯cc . t C ( t ) -2 -4 -6 -8 -10 -12 i (axial-vector) _c c t M e ff a i (axial-vector) _c c FIG. 15: The time-correlation function and the effective mass of the meson interpolator ¯c γ γ i c . C ( t ) -2 -4 -6 -8 -10 -12 ijk j k (tensor)_c c t M e ff a ijk j k (tensor)_c c FIG. 16: The time-correlation function and the effective mass of the meson interpolator ¯c (cid:15) ijk γ j γ k c . ppendix C: C ( t ) and the effective mass of ¯b Γ c t C ( t ) -5 -10 -15 (pseudoscalar) _b c t M e ff a (pseudoscalar)_b c FIG. 17: The time-correlation function and the effective mass of the meson interpolator ¯b γ c . t C ( t ) -5 -10 -15 i (vector) _b c t M e ff a i (vector)_b c FIG. 18: The time-correlation function and the effective mass of the meson interpolator ¯b γ i c . C ( t ) -5 -10 -15 (scalar) _b c t M e ff a (scalar)_b c FIG. 19: The time-correlation function and the effective mass of the meson interpolator ¯bc . t C ( t ) -5 -10 -15 i (axial-vector) _b c t M e ff a i (axial-vector) _b c FIG. 20: The time-correlation function and the effective mass of the meson interpolator ¯b γ γ i c . C ( t ) -5 -10 -15 ijk j k (tensor)_b c t M e ff a ijk j k (tensor)_b c FIG. 21: The time-correlation function and the effective mass of the meson interpolator ¯b (cid:15) ijk γ j γ k c . ppendix D: C ( t ) and the effective mass of ¯b Γ s t C ( t ) -5 -10 -15 (pseudoscalar) _b s t M e ff a (pseudoscalar)_b s FIG. 22: The time-correlation function and the effective mass of the meson interpolator ¯b γ s . t C ( t ) -5 -10 -15 i (vector) _b s t M e ff a i (vector)_b s FIG. 23: The time-correlation function and the effective mass of the meson interpolator ¯b γ i s . C ( t ) -5 -10 -15 (scalar) _b s t M e ff a (scalar)_b s FIG. 24: The time-correlation function and the effective mass of the meson interpolator ¯bs . t C ( t ) -5 -10 -15 i (axial-vector) _b s t M e ff a i (axial-vector) _b s FIG. 25: The time-correlation function and the effective mass of the meson interpolator ¯b γ γ i s . C ( t ) -5 -10 -15 ijk j k (tensor)_b s t M e ff a ijk j k (tensor)_b s FIG. 26: The time-correlation function and the effective mass of the meson interpolator ¯b (cid:15) ijk γ j γ k s ..