Beauty mesons in $N_f=2+1+1+1 $ lattice QCD with exact chiral symmetry

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, with the pion mass $ \sim 700 $ MeV. The gauge ensemble is generated by hybrid Monte Carlo simulation with the Wilson gauge action for the gluons, and the optimal domain-wall fermion action for the quarks. Using point-to-point quark propagators, we measure the time-correlation functions of quark-antiquark meson interpolators with quark contents $\bar b b$, $\bar b c$, $\bar b s$, and $ \bar c c$, and obtain the masses of the low-lying mesons. They are in good agreement with the experimental values, plus some predictions which have not been observed in experiments. Moreover, we also determine the masses of $(b, c, s)$ quarks.


I. INTRODUCTION
In 2007, we performed the first study of treating valence (u, d, s, c, b) quarks as Dirac fermions in quenched lattice QCD with exact chiral symmetry [1,2]. The low-lying mass spectra of mesons with quark contentsbb,bc,bs, andcc were determined, together with the pseudoscalar decay constants. Some of our results (e.g., the masses of η b and h b ) were theoretical predictions at the time of publication, which turn out to be in good agreement with later experimental results. This asserts that it is feasible to treat (u, d, s, c, b) valence quarks 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 lattice QCD with exact chiral symmetry. This motivates the present study. Since the b quark is heavy, with mass m b ∼ 4500 MeV/c 2 , it requires a fine lattice spacing such that the condition m b a < 1 is well satisfied in order to keep the discretization error under control. On the other hand, to keep the finite-volume error of the light hadrons under control, the lattice size L has to be sufficiently large such that M π L 1. These two constraints (a ∼ 0.033 fm and M π L ∼ 4 − 6) together give the lattice size ∼ 170 4 − 260 4 (see Fig. 1), which is beyond the capability of the present generation of supercomputers. Nevertheless, even before the next generation of Exaflop supercomputers will be available ∼ 2025, one may use a smaller lattice to investigate whether the (b, c, s) quarks with physical masses can be dynamically simulated on the lattice, while keeping u and d quarks heavier than their physical masses.
If the pion mass is kept at ∼ 700 MeV/c 2 , then both constraints M π L > 4 and m b a < 1 can be satisfied by the 40 3 × 64 lattice. For domain-wall fermion with the extent N s = 16 in the fifth dimension, the entire hybrid Monte Carlo (HMC) simulation [3] on the 40 3 × 64 × 16 lattice can be performed by one GPU with at least 19 GB device memory, provided that the exact one-flavor pseudofermion action (EOFA) for domain-wall fermion [4] is used. In this study, we use two Nvidia GTX-TITAN-X GPU cards (each of 12 GB device memory) for each stream of HMC simulation, with the peer-to-peer communication between 2 GPUs through the PCIe bus on the motherboard.
The outline of this paper is as follows. In section 2, we recall the basics of lattice QCD with exact chiral symmetry, and discuss what is a viable framework to perform HMC simulation of lattice QCD with both heavy and light domain-wall quarks such that all topological sectors are sampled ergodically and also the chiral symmetry can be peserved to a high precision, i.e., the residual mass of any heavy/light quark flavor is negligible in comparison with its bare mass. In section 3, we describe our lattice setup. In section 4, we determine the low-lying mass spectra of mesons with valence quark contentsbb,bc,bs andcc. In section 4, we determine the masses of (b, c, s) quarks. In section 5, we conclude with some remarks.

A. Preliminaries
Since all quarks in QCD are excitations of Dirac fermion fields, it is vital to preserve this essential feature in lattice QCD. The most theoretically appealing lattice fermion scheme is the domain-wall/overlap fermion [5][6][7], which preserves the exact chiral symmetry at finite lattice 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-wall fermion [8], of which the lattice fermion operator can be written as [D(m q )] xx ;ss (m q ) = (ω s D w + 1) xx δ ss + (ω s D w − 1) xx L ss , where {ω s , s = 1, · · · , N s } are the exact solutions such that the effective 4-dimensional lattice Dirac operator possesses the optimal chiral symmetry for any finite N s . The indices x and x denote the lattice sites on the 4-dimensional lattice, and s and s the indices in the fifth dimension, while the Dirac and color indices have been suppressed. Here D w is the standard Wilson Dirac operator plus a negative parameter −m 0 (0 < m 0 < 2) (m 0 is usually called the domain-wall height), where U µ (x) denotes the link variable pointing from x to x+μ. The operator L is independent of the gauge field, and it can be written as and where m q is the bare quark mass, and m P V = 2m 0 is the Pauli-Villars mass for the optimal DWF. Note that the matrices L ± satisfy L T ± = L ∓ , and R 5 L ± R 5 = L ∓ , where R 5 is the reflection operator in the fifth dimension, with elements (R 5 ) ss = δ s ,Ns+1−s . Thus R 5 L ± is real and symmetric.
Then the pseudofermion action for the optimal DWF can be written as 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 function gives the fermion determinant of the effective 4-dimensional lattice Dirac operator D Ns (m q ), i.e., In the limit N s → ∞, S Ns (H w ) → H w / H 2 w , and D Ns (m q ) goes to In the massless limit m q = 0, D(0) is equal to the overlap-Dirac operator [6], and it satisfies the Ginsparg-Wilson relation [9] where the chiral symmetry is broken by a contact term, i.e., the exact chiral symmetry at finite lattice spacing. Note that (3) does not guarantee that any Ginsparg-Wilson Dirac operator 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 the topological 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] satisfies the 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 operator is the only lattice Dirac operator to possesss topologically exact zero modes satisfying the Atiyah-Singer index theorem on a finite lattice.
However, to perform HMC simulation of lattice QCD with the overlap Dirac operator is prohibitively expensive even for a small lattice (e.g., 16 3 × 32), since it requires to compute the change of the number of exact zero modes n ± at each step of the molecular dynamics [12].
Moreover, the discontinuity of the fermion determinant at the topological boundary highly suppresses the crossing rate between different topological sectors, thus renders HMC failing to sample all topological sectors ergodically. These difficulties can be circumvented by using DWF with finite N s . Firstly, any positive lattice Dirac operator satisfying γ 5 -hermiticity (γ 5 Dγ 5 = D † ) possesses a positive-definite pseudofermion action, without explicit dependence on n ± . Secondly, the step function of the fermion determinant at the topological boundary can be smoothed out by using DWF with finite N s (e.g., N s = 16), then the HMC on the 5-dimensional lattice can sample all topological sectors ergodically and also keep the chiral symmetry to a high precision with the optimal DWF [8,13]. This has been demonstrated for N f = 2 [14], N f = 1 + 1 [4], N f = 2 + 1 + 1 [15], and also N f = 2 + 1 + 1 lattice QCD at the physical point [16].

B. Domain-wall fermion for heavy and light quarks
In this subsection, we discuss which variant of DWF is more capable in capturing the quantum 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 m P V . Otherwise, if m q ∼ m P V , then det(m q )/ det(m P V ) ∼ 1, the internal quark loops are highly suppressed, and the quantum fluctuations of the quark field become mostly quenched. In general, the Pauli-Villars mass is is more severe than the common constraint m q a < 1 for all lattice fermions. In other words, the Shamir/Möbius DWF is not well-suited for studying lattice QCD with heavy quarks.
On the other hand, for the Borici [19]/Optimal [8] DWF, d = 0 and m P V a = 2m 0 1, thus provides the highest ceiling for accommodating the heavy quarks on the lattice, as well as the minimal lattice artifacts due to the mass cutoff. This can be seen by comparing the eigenvalues of their effective 4D Dirac operators in the limit N s → ∞, which is exactly equal to the overlap Dirac operator with the kernel H = cH w (1 + dγ 5 H w ) −1 in the sign function, where c = d = 1/2 for the Shamir/Möbius DWF, while c = 1 and d = 0 for the Borici/Optimal DWF. The eigenvalues of (4) are lying on a circle in the complex plane with radius (m P V − m q )/2, and center at m q + (m P V − m q )/2 on the real axis. are not restricted to a very small circle even for the heavy quark. Moreover, in the chiral limit (left panel), the radius of the eigenvalue circle for the Borici/Optimal DWF is more than 2 times of that of the Shamir/Möbius DWF. This implies that the Borici/Optimal DWF is more capable than the Shamir/Möbius DWF in capturing the short-distance quantum fluctuations 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 what is the optimal chiral symmetry for any finite N s in the fifth dimension, in the sense how its effective 4D lattice Dirac operator can be exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac operator. The exact solution to this problem is given in Ref. [8], with the optimal {ω s } where sn(v s ; κ ) is the Jacobian elliptic function with argument v s (see Eq. (13) in Ref. [8]) and modulus κ = 1 − λ 2 min /λ 2 max . Then S Ns (H w ) is exactly equal to the Zolotarev optimal rational approximation of H w / H 2 w , i.e., the approximate sign function S Ns (H w ) , with degree (n − 1, n) for N s = 2n. Nevertheless, the optimal weights {ω s } in (5) do not satisfy the R 5 symmetry (ω s = ω Ns−s+1 ) which is required for the exact one-flavor pseudofermion action for DWF [4]. The optimal {ω s } satisfying R 5 symmetry is obtained in Ref. [13]. For N s = 2n, the optimal {ω s } satisfying R 5 symmetry are written as where sn(u; κ ) is the Jacobian elliptic function with modulus κ = 1 − λ 2 min /λ 2 max , and K is the complete elliptic function of the first kind with modulus κ . Then the approximate does not satisfy the criterion that the maxima and 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) are preserved, 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 2d Z (0). This can be regarded as the generalized optimal rational approximation (with a constant shift).

III. GENERATION OF THE GAUGE ENSEMBLE
In this section, we give the details of the actions, the algorithms, and the parameters to perform the HMC simulations in this study. Moreover, for the initial 257 trajectories generated by a single node (with 2 Nvidia GTX-TITAN-X GPU cards), the topological charge fluctuation is measured, and the HMC characteristics are presented. Details of the lattice setup are given as follows.

A. The actions
In the following, we present the details of the fermion actions and the gauge action in our HMC simulations.
As noted in Ref. [15], for domain-wall fermions (DWF), to simulate N f = 2 + 1 + 1 amounts to simulate N f = 2 + 2 + 1. Similarly, to simulate N f = 2 + 1 + 1 + 1 amounts to one with det D(m s )/ det D(m P V ) on the LHS. Thus, we perform the HMC simulation with the expression on the RHS of Eq. (7).
For the two-flavor parts, det D( use the N f = 2 pseudofermion action which has been using since 2011 [14], and it can be written as where and L(m q ) is defined in (1) and (2). Here ω ≡ diag{ω 1 , ω 2 , · · · , ω Ns } is a diagonal matrix in the fifth dimension, and D EO/OE w denotes the part of D w with gauge links pointing from even/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 = 2 pseudofermion action (8) is replaced with 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 the exact one-flavor pseudofermion action (EOFA) for DWF [4]. For the optimal DWF, it can be written as ( where φ ± and φ † ± are pseudofermion fields (each of two spinor components) on the 4dimensional lattice, and Here For the gluon fields, we use the Wilson plaquette gauge action [22] at β = 6/g 2 0 = 6.70.
where g 0 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 bare masses of (b, c, s) are tuned to (m b , m c , m s ) = (0.85, 0.20, 0.15) such that they give the masses of the vector mesons Υ(9460), J/ψ(3097), and φ(1020) respectively. The algorithm for 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-wall fermion has been presented in Refs. [4,21]. In the molecular dynamics, we use the Omelyan integrator [23], the multiple-time scale method [24], and the mass-preconditioning [20]. In the following, we summarize the HMC characteristics of the first 257 trajectories. In  Table I.

C. Topological charge fluctuations
In this subsection, we examine the evolution of the topological charge Q t in the first 257 trajectories, and obtain the histogram of its distribution.
In lattice QCD with exact chiral symmetry, the topological charge Q t can be measured by the index of the massless overlap-Dirac operator, since its index satisfies the Atiyah-Singer 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 spherical region of root-mean-square radius R rms = √ 8t, where t is the flow-time. In this study, the flow equation is numerically integrated from t = 0 with ∆t/a 2 = 0.01, and measure the Q clover at t/a 2 = 0.4 which amounts to averaging the gauge field over a spherical region of root-mean-square radius R rms = √ 8t ∼ 1.8a. Then each gauge configuration becomes very 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 the where B −1 x,s;x ,s = δ x,x (P − δ s,s + P + δ s+1,s ) with periodic boundary conditions in the fifth dimension. Then the solution of (12) gives the valence quark propagator 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 mass according to [31], where (D c + m q ) −1 denotes the valence quark propagator with m q equal to the sea-quark mass, tr denotes the trace running over the color and Dirac indices, and the brackets · · · denote the averaging over the gauge ensemble. In the limit N s → ∞, D c is exactly chiral symmetric 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 not well-defined for each gauge configuration, but its impact on any physical observable can be roughly estimated by the difference due to changing the valence quark mass from m q to m q + m res . 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 than 0.007 MeV, which should be negligible in comparison with other systematic uncertainties.
For the vector meson, we average over i = 1, 2, 3 components, namely, 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-correlation function.C The time-correlation function (TCF) and the effective mass of the meson interpolators bΓb,cΓc,bΓc, andbΓ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 which have been measured precisely in high energy experiments.
Our results of the mass spectrum of the low-lying states of bottomonium are summarized in Table III. The time-correlation function and the effective mass ofbΓb are plotted in Appendix A.
The first column in Table III is the Dirac matrix used for computing the time-correlation function (15). The second column is J P C of the state. The third column is the [t 1 , t 2 ] used for fitting the data of C Γ (t) to the usual formula 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 the meson state, where the first error is statistical, and the second is systematic. Here the statistical error is estimated using the jackknife method with the bin size of which the the statistical error saturates, while the systematic error is estimated based on all fittings satisfying χ 2 /dof < 1.2 and |t 2 − t 1 | ≥ 6 with t 1 ≥ 10 and t 2 ≤ 32. The last column is the experimental state we have identified, and its PDG mass value [32].
The analysis and the descriptions in the above paragraph apply to all results obtained in this work, as given in Table III-VI.   TABLE III: The masses of low-lying bottomonium states obtained in this work. The fifth column is 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 a detailed description of each column, see the paragraph with Eq. (16). Evidently, the masses of bottomonium in Table III are   Next, we turn to the charmonium states extracted from the gound states ofcΓc. Our results of the masses of the low-lying states of charmonium are summarized in Table IV. The time-correlation function and the effective mass ofcΓc are plotted in Appendix B.
Evidently, the theoretical masses of charmonium in Table IV are in good agreement with the PDG values. Note that the theoretical result of the hyperfine splitting (1 3 S 1 − 1 1 S 0 ) is 123(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 Table   V. The time-correlation function and the effective mass ofbΓs are plotted in Appendix D. Here we have identified the scalarbs meson with the state B * sJ (5850) observed in high energy 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 energy experiments, thus it serves as a prediction of N f = 2 + 1 + 1 + 1 lattice QCD. Finally, we turn to the heavy mesons with beauty and charm. In Table VI,    In general, Z m should be determined nonperturbatively. However, in this study, the lattice spacing is rather small (a 0.03 fm), thus it is justified to use the one-loop perturbation formula [33] Z s (µ) = 1 + g 2 4π 2 ln(a 2 µ 2 ) + 0.17154 (m 0 = 1.30).
Finally we turn to the strange quark mass. Using (17), the strange quark bare mass m s = 0.0150(2)a −1 is transcribed to where the error bar combines (in quadrature) the statistical and the systematic ones from the lattice spacing and the s quark bare mass. Our result of the strange quark mass (20) is slightly smaller than the PDG value (92.9 ± 0.7) MeV for N f = 2 + 1 + 1 lattice QCD [32].

VI. CONCLUDING REMARK
This study demonstrates that the Dirac b quark can be simulated dynamically in lattice QCD, together with the (c, s, d, u) quarks. Even with unphysically heavy u and d quarks in the sea, the low-lying mass spectra of mesons with valence quark contentsbb,bc,bs, and cc are in good agreement with the experimental values. Also, we have several predictions which 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 and the J P of the pseudovector B s meson state (see Table V). Moreover, we have determined the 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 (∼ 200 4 ) lattice, with the Exaflops supercomputers which will be available ∼ 2025. Then physical observables with any (u, d, s, c, b) quark contents can be computed from the first principles of QCD. This will provide a viable way to systematically reduce the uncertainties in the theoretical predictions of the Standard Model (SM), which are largely stemming from the sector of the strong interaction [34]. This is crucial for unveiling any new physics beyond the standard model (SM), by identifying any discrepancies between the high energy experimental results and the theoretical values derived from the first principles of the SM with all quarks (heavy and light) as Dirac fermions, without using non-relativistic approximation or heavy quark effective field theory for b and c quarks.