Mass Spectra of $D_s$ and $\Omega_c$ in Lattice QCD with $N_f= 2+1+1$ Domain-Wall Quarks

We perform hybrid Monte Carlo simulation of lattice QCD with $N_f=2+1+1 $ optimal domain-wall quarks on the $32^3 \times 64 $ lattice with lattice spacing $a \sim 0.06$ fm, and generate a gauge ensemble with physical $s$ and $c$ quarks, and pion mass $\sim 280 $ MeV. Using 2-quark (meson) and 3-quark (baryon) interpolating operators, the mass spectra of the lowest-lying states containing $s$ and $c$ quarks ($D_s$ and $\Omega_c$) are extracted \cite{Chen:2017kxr}, which turn out in good agreement with the high energy experimental values, together with the predictions of the charmed baryons which have not been observed in experiments. For the five new narrow $\Omega_c$ states observed by the LHCb Collaboration \cite{Aaij:2017nav}, the lowest-lying $\Omega_c(3000)$ agrees with our predicted mass $3015(29)(34)$ MeV of the lowest-lying $\Omega_c$ with $J^P = 1/2^{-}$. This implies that the $ J^P $ of $ \Omega_c(3000) $ is $ 1/2^- $.


Introduction
One of the main objectives of lattice QCD is to extract the hadron mass spectra from the first principles of QCD nonperturbatively. To this end, the hadron mass spectra have to be obtained in a framework which preserves all essential features of QCD, i.e., lattice QCD with overlap/domain-wall fermion, 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 the new physics, or just the approximations (e.g., HQET, NRQCD, partially quenched approximation, etc.) one has used.
In Ref. [1], we use a GPU cluster of 64 Nvidia GTX-TITAN GPUs, and perform the first dynamical simulation of lattice QCD with N f = 2 + 1 + 1 domain-wall quarks on the 32 3 × 64 lattice with extent N s = 16 in the fifth dimension, with physical s and c quarks. To accommodate the c quark without large discretization error, we use a fine lattice (with a ∼ 0.063 fm) such that m c a = 0.55 < 1. Also, to avoid large finite-volume error, we choose the pion mass M π ∼ 280 MeV such that M π L > 3. Even with unphysically u/d quarks in the sea, the mass spectra of hadrons containing c and s quarks turn out in good agreement with experimental results.
In this talk, I review the mass spectra of D s mesons and Ω c baryons obtained in Ref. [1], and discuss the physical implications, and also predictions in high energy experiments.

Lattice Setup
As pointed out in Ref. [1], for the domain-wall fermion, to simulate N f = 2 + 1 + 1 amounts to simulate N f = 2 + 2 + 1, since where D(m q ) denotes the domain-wall fermion operator with bare quark mass m q , and m PV the mass of the Pauli-Villars field. Since the simulation of 2-flavors is most likely faster than the simulation of one-flavor, it is better to simulate N f = 2 + 2 + 1 than N f = 2 + 1 + 1.

Action and simulation
For the gluon fields, we use the Wilson plaquette gauge action at β = 6/g 2 0 = 6.20. For the quark fields, we use the optimal domain-wall fermion actions [3,4]. For the HMC simulation of the 2-flavors, we use the pseudofermion action for 2-flavors lattice QCD with domain-wall fermion as defined by Eq. (14) in Ref. [5]. For the simulation of the one-flavor, we use the exact one-flavor pseudofermion action (EOFA) for domain-wall fermion, as defined by Eq. (23) in Ref. [6]. The parameters of the pseudofermion actions are fixed as follows. For the domain-wall fermion operator D(m q ) defined in Eq. (2) of Ref. [5], we fix c = 1, d = 0 (i.e., H = H w ), m 0 = 1.3, N s = 16, and λ max /λ min = 6.20/0.05. For the 2-flavors action, the optimal weights {ω s , s = 1, · · · , N s } are computed according to Eq. (12) in Ref. [3] 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. [4], which are the optimal weights satisfying the R 5 symmetry, giving the approximate sign function S (H) of the effective 4D Dirac operator satisfying the bound 0 . 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 simulation of 2-flavors has been outlined in Ref. [5], while that for the exact one-flavor action (EOFA) has been presented in Ref. [6].
We generate the initial 460 trajectories with two Nvidia GTX-TITAN cards. After discarding the initial 300 trajectories for thermalization, we sample one configuration every 5 trajectories, resulting 32 "seed" configurations. Then we use these seed configurations as the initial configurations for 32 independent simulations on 32 nodes, each of two Nvidia GTX-TITAN cards. Each node generates 50 − 85 trajectories independently, and all 32 nodes accumulate a total of ∼ 2000 trajectories. From the saturation of the binning error of the plaquette, as well as the evolution of the topological charge, we estimate the autocorrelation time to be around 5 trajectories. Thus we sample one configuration every 5 trajectories, and obtain a total of 400 configurations for physical measurements.
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 [10] 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 (1) gives the valence quark propagator Each column of the quark propagator is computed with 2 Nvidia GTX-TITAN GPUs in one computing node, attaining more than one Teraflops/sec (sustained).

Residual masses
To measure the chiral symmetry breaking due to finite N s , we compute the residual mass according to Eq. (45) in Ref. [11]. 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 1. 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.

Mass spectra of D s mesons and Ω c baryons
We construct quark-antiquark interpolators for mesons, and 3-quark interpolators for baryons, and measure their time-correlation functions using the point-to-point quark propagators computed with the same parameters of the sea quarks. Then we extract the mass of the lowest-lying hadron states from the time-correlation function, following the procedures outlined in Refs. [12][13][14].

Mass spectrum of D s mesons
The time-correlation function of the D s meson interpolatorcΓs is measured according to the formula for scalar (S ), pseudoscalar (P), vector (V), axial-vector (A), and tensor (T ), with Dirac matrix Γ = {1 I, γ 5 , γ i , γ 5 γ i , γ 5 γ 4 γ i = i jk γ j γ k /2} respectively. Note that the Dirac bilinear covariantq i jk γ j γ k q is  often called as "tensor" in the textbook. However, it transforms like axial-vector since its J P = 1 + , different from the usual terminology "tensor meson" which refers to the mesons with J = 2. In the following "tensor meson" always refers to that with Γ = i jk γ j γ k and J P = 1 + . For the vector meson, we average over i = 1, 2, 3 components. Similarly, we perform the same averaging for the axial-vector and the tensor mesons. Moreover, to enhance statistics, we average the forward and the backward time-correlation function.
The time-correlation functions of all meson channels are plotted in the left panel of Fig. 1. The effective mass of the scalar (Γ = 1) is plotted in the right panel of Fig. 1, and those of the axial vector (Γ = γ 5 γ i ) and tensor (Γ = i jk γ j γ k ) are plotted in Fig. 2. Since bothcγ 5 γ i s andc i jk γ j γ k s have J P = 1 + , one expects that there are mixings between them. However, from the time-correlation functions (see the left panel of Fig. 1), they seem to be two distinct states with different masses, with little overlap. Thus we extract the masses of the lowest-lying states of thecs mesons from each channel (Γ) individually, and the results are summarized in Table 2. The first column is the Dirac matrix. The second column is J P of the state. The third column is the time interval [t 1 , t 2 ] for fitting the data of the time-correlation function C Γ (t) to the formula to extract the meson mass M and the amplitude z = | H|QΓq|0 |, where H denotes the lowest-lying meson state with zero momentum, and the excited states have been neglected in (2). The fifth column is the mass M of the state, where the first error is statistical, and the second is systematic error. Here the statistical error is estimated using the jackknife method with the bin-size of which the statistical error saturates, while the systematic error is estimated based on all fittings satisfying χ 2 /dof ≤ 1.1 and |t 2 − t 1 | ≥ 5 with t 1 ≥ 10 and t 2 ≤ 30. The last column is the corresponding state in high energy experiments, with the PDG mass value [15]. Evidently, our results of the mass spectrum of the lowestlying states of the the D s mesons are in good agreement with the PDG values. This implies 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. 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, the time-correlation function of any D s quark-antiquark interpolator is well fitted to the form of single meson state (2) on a plateau with |t 1 − t 2 | ≥ 5. This implies that are much less than one for t ∈ [10,54], and the nearby scattering states have little overlap with the physical meson. Following our notations in Ref. [12], the interpolating operators for Ω c baryons are [c(Cγ 5 )s]s and (cCγ µ s)s, where C is the charge conjugation operator. The time-correlation function of any baryon interpolator B is defined as C αβ (t) = x B xαB0β , which can be expressed in terms of quark propagators. The ensemble-averaged time-correlation function is fitted to the usual formula

Mass spectrum of Ω
where m ± are the masses of even and odd parity states. Thus, one can use the parity projector P ± = (1 ± γ 4 )/2 to project out two amplitudes, For sufficiently large T , there exists a range of t such that, in A ± , the contributions due to the opposite parity state are negligible. Thus m ± and Z ± can be extracted by a single exponential fit to A ± = Z ± e −m ± at , for the range of t in which the effective mass m eff (t) = ln(A ± (t)/A ± (t + 1)) attains a plateau. 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 k j (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 i j (t). To enhance the statistics, we use 3 i=1 C 3/2 ii (t)/3 to extract the mass of the J = 3/2 state.
In Table 3, we summarize the masses of Ω and Ω c baryon states obtained in Ref. [1]. The mass value in the fifth column is obtained by correlated fit, where the first error is statistical, and the second is systematic error. Here the statistical error is estimated using the jackknife method with the bin-size of which the statistical error saturates, while the systematic error is estimated based on all fittings satisfying χ 2 /dof ≤ 1.2 and |t 2 − t 1 | ≥ 5 with t 1 ≥ 10 and t 2 ≤ 30. Evidently, the masses of Ω(3/2 + ), Ω(3/2 − ), Ω c (1/2 + ), and Ω c (3/2 + ) are in good agreement with the PDG values in the last column. For Ω c (1/2 − ) and Ω c (3/2 − ), they had not been observed in experiments when Ref. [1] was published in January 2017. In March 2017, five new narrow Ω c states were observed by the LHCb Collaboration [2], the lowest-lying Ω c (3000) agrees with our predicted mass 3015(29)(15) MeV of the lowest-lying Ω c with J P = 1/2 − . This implies that the J P of Ω c (3000) is 1/2 − .

Summary and Outlook
In Ref. [1], 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 (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 2-3. 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. Currently, we are generating additional 2 gauge ensembles with pion masses ∼ 320 − 400 MeV, which can be used for extrapolation to the physical limit. 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 much less than our estimated statistical and systematic errors. For thecs meson states in Table 2, 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 Ω and Ω c in Table 3, the masses of Ω(3/2 + ), Ω(3/2 − ), Ω c (1/2 + ), and Ω c (3/2 + ) are in good agreement with the PDG values. For Ω c (1/2 − ) and Ω c (3/2 − ), they had not been observed in experiments when Ref. [1] was published in January 2017. In March 2017, five new narrow Ω c states were observed by the LHCb Collaboration [2], the lowest-lying Ω c (3000) agrees with our predicted mass 3015(29)(34) MeV of the lowest-lying Ω c with J P = 1/2 − . This implies that the J P of Ω c (3000) is 1/2 − . Now the challenge is to find out the full spectrum of Ω c (including the excited states) in the framework of lattice QCD with domain-wall quarks, and to see whether they can be identified with the five new narrow Ω c states observed by the LHCb Collaboration.