Lattice QCD with $N_f = 2+1+1 $ domain-wall quarks

We perform hybrid Monte Carlo simulation of (2+1+1)-flavors lattice QCD with the optimal domain-wall fermion (which has the effective 4D Dirac operator exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac operator). The gauge ensemble is generated on the $32^3 \times 64 $ lattice with the extent $ N_s = 16 $ in the fifth dimension, and with the plaquette gauge action at $ \beta = 6/g^2 = 6.20 $. The lattice spacing ($ a \simeq 0.063 $ fm) is determined by the Wilson flow, using the value $ \sqrt{t_0} = 0.1416(8) $ fm obtained by the MILC Collaboration for the $(2+1+1)$-flavors QCD. The masses of $s$ and $c$ quarks are fixed by the masses of the vector mesons $ \phi(1020) $ and $ J/\psi(3097) $ respectively; while the mass of the $u/d$ quarks is heavier than their physical values, with the unitary pion mass $ M_\pi \simeq 280$ MeV (and $ M_\pi L \simeq 3 $). We compute the point-to-point quark propagators, and measure the time-correlation functions of meson and baryon interpolators. Our results of the mass spectra of the lowest-lying hadrons containing $ s $ and $ c $ quarks are in good agreement with the high energy experimental values, together with the predictions of the charmed baryons which have not been observed in experiments.


I. INTRODUCTION
Since the discovery of the Higgs scalar in 2012, the Standard Model (SM) emerged in mid 1970s looks to be complete in the sense that all major predictions of the SM have been realized in high energy experiments, and almost all high energy experimental data can be understood in the framework of SM, except the matter-antimatter asymmetry and the origin of the neutrino masses. Currently, the challenge of high energy physics is to find out whether there is any new physics beyond the SM, in view of the generation puzzle and the large number of parameters in the SM, which suggest that the SM is an effective theory at the scale probed by the present generation of high energy accelerators. In order to identify any discrepancies between the high energy experimental results and theoretical values derived from the SM, the latter have to be obtained in a framework which preserves all essential features of the SM. Otherwise, it is difficult to determine whether such a discrepancy is due to new physics, or just the approximations (or models) one has used. So far, the largest uncertainties in the theoretical predictions of the SM stem from the sector of the strong interaction, namely, QCD. Theoretically, lattice QCD is the most viable framework to tackle QCD nonperturbatively from the first principles. However, in practice, it is difficult to simulate dynamical u, d, s, c, and b quarks at their physical masses (ranging from ∼ 3 − 4500 MeV), in a sufficiently large volume and small enough lattice spacing such that the finite-volume and discretization errors are both well under control. Note that the t quark can be neglected in QCD simulations since it is extremely short-lived and it decays to W-boson and b/s/d quarks before it can interact with other quarks through the gluons.
Even after neglecting the t quark, to simulate u, d, s, c, and b quarks at their physical masses is still a very challenging problem. For example, if one designs the simulation close to the physical pion mass with M π ≃ 140 MeV and M π L > 6 (to keep finite-volume error under control), then it would require a lattice of size ∼ 100 4 to accommodate physical c quark with sufficiently small discretization error, not to mention the much heavier b quark. The current generation of supercomputers with ∼ 100 Petaflops seems to be marginal for this purpose, and the next generation of supercomputers with Exaflops is required to simulate (u, d, s, c, b) quarks at their physical masses.
With our rather limited resources, we can only afford to perform lattice QCD simulations with domain-wall quarks on a 32 3 × 64 lattice, using a GPU cluster with 64 Nvidia GTX-TITAN GPUs. Now, even after neglecting the dynamical b quark, we still have two options.
One way is to neglect the c quark, and simulate (2 + 1)-flavors QCD on a coarse lattice such that M π ∼ 140 MeV and M π L > 3, for studying the phenomenology involving the light quarks. Instead, we simulate (u, d, s, c) quarks with sufficiently fine lattice spacing satisfying m c a < 1, for studying the charm physics, which in turn must render the unitary pion [m u/d (valence) = m u/d (sea)] heavier than 140 MeV such that M π L > 3 to avoid large finite-volume error.
Even with unphysically heavy u/d quarks in the sea, the mass spectra of hadrons containing c and s quarks may turn out to be in good agreement with high energy experimental results, as we have observed in our previous studies, for N f = 0 lattice QCD [1,2] and In this paper, we examine to what extent this scenario holds for (2 + 1 + 1)-flavors QCD.
We perform hybrid Monte Carlo simulation of lattice QCD with N f = 2 + 1 + 1 optimal domain-wall quarks [4,5]  First, we point out that, for the domain-wall fermion, to simulate N f = 2 + 1 + 1 amounts to simulate N f = 2 + 2 + 1, according to the identity where D(m q ) denotes the domain-wall fermion operator with bare quark mass m q , and m P V the mass of the Pauli-Villars field. Since the simulation of 2-flavors is much faster than the simulation of one-flavor, it is better to simulate N f = 2 + 2 + 1 than N f = 2 + 1 + 1.
Whether (2) is more efficient than (3) or vice versa depends on the computational platform, the algorithm, the lattice size, and the parameters of the action. In this work, we choose (3) for our HMC simulations.
For the gluon fields, we use the Wilson plaquette gauge action at β = 6/g 2 0 = 6.20. For the two-flavor parts, we use the pseudofermion action for 2-flavors lattice QCD with the optimal domain-wall quarks, as defined in Eq. (14) of Ref. [6]. For the one-flavor part, we use the exact pseudofermion action for one-flavor domain-wall fermion, as defined by Eq. For the 2-flavors action, ω s are computed according to Eq. (12) in Ref. [4] such that the effective 4D Dirac operator is exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac operator with bare quark mass m q . For the one-flavor action, ω s are computed according to Eq. (9) in Ref. [5], which are the optimal weights satisfying the R 5 symmetry, giving the approximate sign function S(H) of the effective 4D Dirac operator . We perform the HMC simulation of (2+1+1)-flavors QCD on the L 3 × T = 32 3 × 64 lattice, with the u/d quark mass m u/d a = 0.005, the strange quark mass m s a = 0.04, and the charm quark mass m c a = 0.55, where the masses of s and c quarks are fixed by the masses of the vector mesons φ(1020) and J/ψ(3097) respectively. The algorithm for simulating 2-flavors of optimal domain-wall quarks has been outlined in Ref. [6], while the exact one-flavor algorithm (EOFA) for domain-wall fermions has been presented in Ref. [7].
Here we note that EOFA outperforms the rational hybrid Monte Carlo algorithm (RHMC) [8], no matter in terms of the memory consumption or the speed [9,10].
In the molecular dynamics, we use the Omelyan integrator [11], and the Sexton-Weingarten multiple-time scale method [12]. Moreover, we introduce an auxiliary heavy fermion field with mass m H (m q ≪ m H ≪ m P V ) similar to the case of the Wilson fermion [13], the so-called mass preconditioning. For the 2-flavors parts, mass preconditioning is only applied to the u/d quark factor of (3) with m H a = 0.1. For the one-flavor part, a novel mass preconditioning has been devised for the EOFA [14], which is ∼ 20% faster than the mass preconditioning we have used in Refs. [7,9]. Also, based on the fact that in EOFA the fermion force of the φ 1 field is much smaller than that of the φ 2 field, the gauge momentum updating by these two forces can be set at two different time scales. Furthermore, we have developed a generalized multiple-time scale method with the flexibility of assigning an arbitrary updating time interval to any fermion force provided that the updating time interval of the gauge force is an integer multiple of the lowest common multiplier (LCM) of the updating intervals of all fermion forces. This feature is essential for tuning the parameters to attain optimal efficiency. The details of our simulations will be presented in a forthcoming long paper [15].
We generate the initial 460 trajectories with two Nvidia GTX-TITAN cards (each with device memory ≥ 6 GB). After discarding the initial 300 trajectories for thermalization, we sample one configuration every 5 trajectories, resulting 32 "seed" configurations.
, where the matrix-valued field tensor F µν (x) is obtained from the four plaquettes surrounding x on the (μ,ν) plane. Even though this topological charge is not exactly equal to an integer, it gives a qualitative picture to demonstrate that our HMC simulation samples all topological sectors ergodically. For a rigorous determination of the topological charge and susceptibility, it requires to project the zero modes of the overlap Dirac operator [16,17] for each gauge configuration, which is beyond the scope of this paper.
We compute the valence quark propagator of the 4D effective Dirac operator with the point source at the origin, and with the mass and other parameters exactly the same as those of the sea quarks. First, we solve the following linear system with mixed-precision conjugate gradient algorithm, for the even-odd preconditioned D (see Eq. (12) in Ref. [6]), 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 (5) gives the valence quark propagator Each column of the quark propagators is computed with 2 Nvidia GTX-TITAN GPUs in one computing node, attaining more than one Teraflops/sec (sustained).  To measure the chiral symmetry breaking due to finite N s , we compute the residual mass according to the formula [21] m res = tr(D c + m q ) −1 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 · · · U denote the averaging over the gauge ensemble. For the 400 gauge configurations generated by HMC simulation of lattice QCD with N f = 2 + 1 + 1 optimal domain-wall quarks, the residual masses of u/d, s, and c quarks are listed in Table I. We see that the residual mass of the u/d quark is ∼ 1.2% of its bare mass, amounting to 0.19(4) MeV, which is expected to be much smaller than other systematic uncertainties. The residual masses of s and c quarks are even smaller, 0.11(3) MeV, and 0.07(3) MeV respectively.

III. MASS SPECTRUM OF HADRONS CONTAINING s AND c QUARKS
One of the main objectives of lattice QCD is to extract the mass spectrum of QCD nonperturbatively from the first principles. Even though our N f = 2 + 1 + 1 gauge ensemble is generated with unphysically heavy u/d quarks (with M π ≃ 280 MeV), we suspect that the mass spectrum and the decay constants of hadrons containing c and s quarks may turn out to be in good agreement with high energy experimental results, as we have observed in N f = 0 lattice QCD [1,2] and N f = 2 lattice QCD [3] respectively. In the following, we examine to what extent this scenario is realized in the spectrum of lattice QCD with N f = 2 + 1 + 1 optimal domain-wall quarks.
Following our previous studies [1][2][3], we construct quark-antiquark interpolators for mesons, and 3-quarks interpolators for baryons, and measure their time-correlation functions    For thecs meson states in Table III,  quarks in the sea. It is interesting to see that the masses of the scalar meson D * s 0 (2317), and the axial-vector mesons D s1 (2460) and D s1 (2536) can be obtained with quark-antiquark interpolators, without invoking 4-quark interpolators like DK and D * K. We note that a recent study [23] of N f = 2 + 1 lattice QCD with nonperturbatively improved Wilson-clover fermions and the same fermion action for the valence quarks, using quark-antiquark interpolators, also obtained the masses of the lowest-lyingcs meson states compatible with the experimental values.
Note that in the physical limit, D * s 0 (2317) is about 41 MeV below the DK threshold, and D s1 (2460) is 44 MeV below the D * K threshold, while D s1 (2536) is 32 MeV above the D * K threshold. Thus it seems to be necessary to consider the effects of the nearby scattering states, e.g., by incorporating 4-quark interpolators like DK and D * K. However, for our gauge ensemble, the DK threshold is about 156 MeV above thecs scalar meson state, and the D * K threshold is more than 220 MeV and 146 MeV above thecs axial-vector meson states. Moreover, since the time-correlation function is well fitted to the form of single meson state (6) on plateaus with |t 1 − t 2 | ≥ 6, this implies that the ratios are much less than one (at least for t ∈ [10, 54]), for our gauge ensemble. Nevertheless, it is still interesting to check whether the masses of these states would be affected by the threshold effects, by incorporating 4-quark interpolators DK and D * K, and performing variational analysis on the correlation matrices of both 2-quark and 4-quark interpolators, similar to the study in Ref. [24], especially for the gauge ensembles approaching the physical limit.
Next, we turn to the baryons with s and c quarks, Ω, Ω c , Ω cc , and Ω ccc . Following the notations in our previous study [1], For baryon interpolating operator like B µ = (q 1 Cγ µ q 2 )q 3 , spin projection is required to extract the J = 3/2 state, since it also overlaps with the J = 1/2 state. The spin J = 3/2 projection for the time-correlation function reads where C kj (t) = x B k ( x, t)B j ( 0, 0) . Then the mass of the J = 3/2 ± state can be extracted from any one of the 9 possibilities (i, j = 1, 2, 3) of C 3/2 ij (t). To enhance the statistics, we ii (t)/3 to extract the mass of the J = 3/2 state. Following the procedures outlined in our previous study [1], we obtain the masses of Ω, Ω c , Ω cc and Ω ccc , as summarized in Table IV. The mass value in the fifth column is obtained by correlated fit, where the first error is statistical, and the second is systematic error.
It is interesting to point out that the mass of Ω c (3/2 + ) was predicted to be 2756(32) MeV in our quenched study [1], before it was observed by the Belle Collaboration in 2009, with the measured mass 2765.9 ± 2.0 MeV [25]. In other words, for lattice QCD with exact chiral symmetry, even in the quenched approximation, it can give the mass spectra of heavy hadrons reliably. This scenario also holds for heavy mesons, e.g., in our quenched study of mesons containing b, c, and s quarks [2], we predicted the mass of η b to be 9383(4)(2) MeV, before η b was discovered by the BaBar Collaboration in 2008 [26], with the measured mass 9388.9 +3.1 −2.3 ± 2.7 MeV. We note that there are several recent lattice studies of the mass spectra of charmed baryons (see, e.g., Refs. [23,[27][28][29][30]), in the framework of N f = 2, 2+1, and 2 + 1 + 1 lattice QCD, with different fermion actions for the c quark, and/or the c quark is absent in the sea.
A detailed review of lattice results of charmed baryons is beyond the scope of this paper.

IV. SUMMARY AND CONCLUDING REMARKS
In this paper, we present the first study of lattice QCD with N f = 2 + 1 + 1 domain-wall quarks. Using 64 Nvidia GTX-TITAN GPUs evenly distributed on 32 nodes, we perform the HMC simulation on the 32 3 × 64 × 16 lattice, with lattice spacing a ∼ 0.063 fm. Even though the mass of u/d quarks is unphysically heavy (with unitary pion mass ∼ 280 MeV), the masses of hadrons containing c and s quarks turn out in good agreement with the experimental values, as summarized in Tables II-IV. However, extrapolation to the physical limit (with M π = 140 MeV) is still required, though we do not expect significant changes in the mass spectra of hadrons containing s and c quarks. Since we have generated only one gauge ensemble, it is impossible for us to perform extrapolation to the physical limit, not to mention taking the continuum limit and the infinite volume limit. Nevertheless, comparing the mass spectra in Tables II-IV to those of N f = 2 lattice QCD on a 24 3 × 48 lattice with a ∼ 0.063 fm [31], we conclude that the finite volume uncertainty is much less than the estimated statistical and systematic errors. About the discretization error, since the lattice spacing (a ∼ 0.063 fm) is sufficiently fine, and our lattice action is free of O(a) lattice artifacts, we expect that the discretization error is also much less than our estimated statistical and systematic errors.
For thecs meson states in Table III, our results show that they are conventional meson states composed of valence quark-antiquark, interacting through the gluons with the quantum fluctuations of (u, d, s, c) quarks in the sea, even for the scalar meson D * s 0 (2317), and the axial-vector mesons D s1 (2460) and D s1 (2536).
For the mass spectra of baryons in Table IV,  To address the challenge of finding out whether there is any new physics beyond the SM, it requires to pin down the theoretical uncertainties largely from the sector of the strong interaction, before one can identify any discrepancies between the experimental results and the theoretical values derived from the SM. To this end, the latter have to be obtained in a framework which preserves all essential features of QCD, i.e., lattice QCD with exact chiral symmetry, and also in the unitary limit (with the valence and the sea quarks having the same masses and the same Dirac fermion action). Otherwise, it is difficult to determine whether any discrepancy between the experimental result and the theoretical value is due to new physics, or just the approximations (e.g., HQET, NRQCD, partially quenched approximation, etc.) one has used.
To conclude, this work asserts that it is feasible to perform large-scale lattice QCD simulations with N f = 2 + 1 + 1 domain-wall quarks, with good chiral symmetry, and sampling all topological sectors ergodically. It provides the ground work for future largescale lattice QCD simulations with dynamical (u, d, s, c, b) domain-wall quarks.