Charmed Tetraquarks Tcc and Tcs from Dynamical Lattice QCD Simulations

Charmed tetraquarks $T_{cc}=(cc\bar{u}\bar{d})$ and $T_{cs}=(cs\bar{u}\bar{d})$ are studied through the S-wave meson-meson interactions, $D$-$D$, $\bar{K}$-$D$, $D$-$D^{*}$ and $\bar{K}$-$D^{*}$, on the basis of the (2+1)-flavor lattice QCD simulations with the pion mass $m_{\pi} \simeq $410, 570 and 700 MeV. For the charm quark, the relativistic heavy quark action is employed to treat its dynamics on the lattice. Using the HAL QCD method, we extract the S-wave potentials in lattice QCD simulations, from which the meson-meson scattering phase shifts are calculated. The phase shifts in the isospin triplet ($I$=1) channels indicate repulsive interactions, while those in the $I=0$ channels suggest attraction, growing as $m_{\pi}$ decreases. This is particularly prominent in the $T_{cc} (J^P=1^+,I=0)$ channel, though neither bound state nor resonance are found in the range $m_{\pi} =410-700$ MeV. We make a qualitative comparison of our results with the phenomenological diquark picture.


Introduction
One of the long standing challenges in hadron physics is to establish and classify genuine multiquark states other than baryons (3 quark states) and mesons (quark-antiquark states) [1]. In particular, charmed tetraquarks (such as T cc (ccūd), T cs (csūd)) and bottomed tetraquarks (such as T bb , T bc , T bs ) are unambiguous candidates for such multiquark states [2,3,4,5], since there are no annihilations among the four quarks. If they form bound states or resonances with respect to the corresponding two-meson thresholds, they could be experimentally observed in B-factories and relativistic heavy-ion colliders [6,7,8,9].
In this Letter, we exclusively investigate the charmed tetraquarks T cc and T cs . To understand a possible reason why they may appear as bound states below the two meson threshold, let us consider the diquark picture [10] as a working hypothesis, whereūd in the color 3, spin-singlet (S=0), isospin-singlet (I=0) channel is denoted as a "good" diquark, due to the large attraction betweenū andd generated through a gluon exchange. We now assume that good diquarks are the main substructure of charmed tetraquarks. Then the color-singletness of hadrons and the quark Pauli principle constrain the possible low-mass tetraquarks as follows.
(i) T cc (J P = 1 + , I = 0) in whichūd forms a good diquark while the diquark cc with color 3 * and S = 1 has a weak repulsion between the two charm quarks. This state couples to the D-D * system. (ii) T cs (J P = 1 + , I = 0) in whichūd forms a good diquark while the diquark cs with color 3 * and S = 1 has a weak repulsion. This state couples to theK-D * system. (iii) T cs (J P = 0 + , I = 0) in whichūd forms a good diquark while the diquark cs with color 3 * and S = 0 has an attraction. This state couples to theK-D system.
Quantitatively, however, predictions for the binding energies of charmed tetraquarks widely spread, ranging from negative values (resonance) to 100 MeV (deeply bound) with respect to the two-meson thresholds, depending on the details of the dynamical models (diquark model, dynamical four-body calculation in the constituent quark model, meson-meson molecular model, lattice QCD in the heavy-quark limit, etc) [11,12,13,14,15,16]. Therefore, a quantitative prediction for charmed tetraquarks requires a careful study in full lattice QCD with a finite charm-quark mass 1 .
In this Letter, we report our first results on the interactions in the D-D,K-D, D-D * andK-D * systems in (2+1)-flavor lattice QCD simulations. The dynamics of the charm quarks are incorporated on the lattice with the relativistic heavy quark action [19]. The meson-meson scattering phase shifts are derived from the corresponding potentials calculated on the lattice by the HAL QCD method [20,21,22] (reviewed in [23]), which was recently shown to be quite accurate to describe some meson-meson scattering phase shifts [24].
This Letter is organized as follows. In Sec. 2, we present the HAL QCD method to extract the potential between two mesons. We then show the numerical setup of our lattice QCD simulations in Sec. 3. In Sec. 4, we show our numerical results for the potentials, from which scattering phase shifts and scattering lengths are extracted for three different quark masses. Sec. 5 is devoted to a summary and a discussion.

HAL QCD method for the meson-meson interaction
In QCD, the two-meson correlation function can be expanded as with A n = n| J h 1 h 2 (t = 0) |0 , J h 1 h 2 (t = 0) stands for a source operator at t src. = 0 which creates two meson states and O h 1,2 is a point-like interpolating sink operator for the hadron h 1,2 . W n = m 2 1 + k 2 n + m 2 2 + k 2 n is the relativistic energy of the n-th eigenstate |n for two mesons, and ellipses represent inelastic contributions.
Consider t sufficiently larger than t src. that the contributions from elastic scattering states and possible bound states remain while those from inelastic states become negligible. Then, Eq. (1) becomes F ( r, t) → n A n φ n ( r)e −Wnt , where φ n ( r) is an equal-time Nambu-Bethe-Salpeter (NBS) wave function [25], from which the HAL QCD potential U is defined [21] as a solution of for all elastic eigenstates n, where H 0 = −∇ 2 /2µ with µ = m 1 m 2 /(m 1 + m 2 ) and E n = k 2 n /2µ is a kinetic energy. Here the non-local but n-independent potential U( r, r ′ ) can be shown to exist, by explicitly constructing it as whereφ * n ( r ′ ) is the dual basis associated with φ n ( r ′ ), and the summation over n is restricted to elastic channels. (For details and proofs, see [21,23]).
In principle, the potential is extracted from F ( r, t) at large t, when it is dominated by the n = 0 state (i.e. the ground state) contribution [21,23]. In practice, however, F ( r, t) is usually noisier at larger t, so that an accurate determination of potentials in this way becomes difficult.
To overcome this practical difficulty, an alternative method has been proposed in [22]. Since U( r, r ′ ) is n-independent by definition, a normalized correlation function R( r, t) = F ( r, t)/e −(m 1 +m 2 )t satisfies where the non-relativistic approximation that In the velocity expansion of the non-local potential, we finally obtain the leading order potential as within the non-relativistic approximation.
To extract S-wave potentials on the lattice, we consider the projection of the normalized correlation function on the A + 1 representation of the cubic group (containing the J = 0 representation of the rotational group) where g ∈ O are elements of the cubic group, and χ (A 1 ) (g)(≡ 1) are the associated characters of the A 1 representation.
As for the charm quark, we employ a relativistic heavy quark (RHQ) action proposed in Ref. [19], which is designed to remove the leading and nextto-leading order cutoff errors associated with heavy quark mass, O((m Q a) n ) and O((m Q a) n (aΛ QCD )), respectively. The RHQ action is given by Parameters of the action are κ Q , r s , r t , ν, c B and c E , while the redundant parameter r t is chosen to be 1. In our simulations, we take the same parameters as in Ref. [29], where the 1S charmonium mass and its relativistic dispersion relation are reproduced. The RHQ parameters are summarized in Table 1.
Periodic boundary conditions are imposed in the three spacial directions, while Dirichlet boundary conditions are taken for the temporal direction at t/a = ±32 to avoid contaminations from the opposite propagation of mesons in time. Throughout this study, we employ local interpolating operators for mesons, φ(x) =q(x)Γq(x), where Γ denotes a 4 × 4 matrix acting on spinor indices. We take Γ = γ 5 for pseudo-scalar mesons (D andK) and Γ = γ i for vector mesons (D * ). Meson masses calculated in this work are listed in Table 2 together with the number of configurations used in this work. We measure the correlation function in Eq. (1) with a source at one time-slice for each configuration, and the forward and backward propagations are averaged to enhance the statistics. We have checked that the dispersion relation of the 1S charmonium state at our heaviest pion mass, m π ∼ 700MeV, gives a reasonable value of the effective speed of light, c eff = 0.987 (2).  As for the source operators of the D-D, D-D * ,K-D andK-D * in isospin I channels, we take the following wall sources:

Meson-meson potentials and scattering phase shifts
With the above setup, we study the S-wave meson-meson interactions in the following channels related to T cc and T cs with J P = 0 + , 1 + : D-D (J P = 0 + , I = 1),K-D (J P = 0 + , I = 0, 1), D-D * (J P = 1 + , I = 0, 1), andK-D * (J P = 1 + , I = 0, 1). We show the potentials calculated from Eq. (5) at the time-slice t/a = 16. Since the energy differences between the elastic and inelastic thresholds are 200-300MeV as shown in Table 3, the inelastic contributions in Eq. (1) are expected to be suppressed in the time-slices t/a > 11.
In Fig. 1, we show our results for the S-wave meson-meson potentials in the I = 1 channels with J P = 0 + (left panels) and J P = 1 + (right panels). We find that all the potentials in the I = 1 channel are repulsive at all distances. This observation is consistent with the absence of goodūd diquarks in the I = 1 channel, as discussed in Section 1. Since the quark mass dependence of the potentials is rather weak, it is unlikely that interactions in these channels turn into strong attractions to form bound states even at the physical quark mass. Phenomenological models also predict the absence of bound states in the I = 1 channels.
In Fig. 2, we show our results for the S-wave meson-meson central potentials in the I = 0 channels with J P = 0 + (left panels) and J P = 1 + (right Threshold energies m π =411(2) m π =572(2) m π =699(1) E DD 3805 (5)    panels). Contrary to the previous results in the I = 1 channels, all the potentials in the I = 0 channels show attractions at all distances without repulsive core. In addition, we find that all potentials become more attractive as the pion mass decreases, and the attraction at short distance (r = 0.2 ∼ 0.3fm) in theK-D channel is stronger than in the D-D * andK-D * channels. Such tendency is again consistent with the existence of goodūd diquarks in the I = 0 channel, as discussed in Section 1.
To investigate the possible existence of bound states or resonances in the I = 0 channels, we fit the potentials with analytic functions of r and solve the Schrödinger equation with the fitted potentials at a given timeslice. We employ a multi-range Gaussian form to fit the potential, namely g(r) ≡ Nmax n=1 V n · exp(−ν n r 2 ) with V n and ν n being fit parameters. For all cases, we obtain good fits with χ 2 /d.o.f. ≃ 0.6 for N max = 4. Repeating this at different time-slices, the mean values of scattering phase shifts are obtained from the weighted average over the time-slices t/a = 13 through 18. Statistical errors for the scattering phase shifts are calculated by the jackknife method, and systematic errors are evaluated by the difference between the weighted average of the phase shifts over the time-slices t/a = 13 − 15 and t/a = 16 − 18. Fig. 3 shows the resultant S-wave scattering phase shifts as a function of the meson-meson center-of-mass energy in the I = 0K-D (J P = 0 + ), K-D * (J P = 1 + ) and D-D * (J P = 1 + ) channels. In Table 4 and Fig. 4, we give the corresponding scattering lengths. We do not find the negativeenergy eigenvalues corresponding to the bound state solutions by solving the Schrödinger equation with the potentials shown in Fig. 2. Fig. 3 also indicates that there are no bound states or resonances in this range of pion masses, m π = 410 ∼ 700 MeV. Although the potentials for D-D * andK-D * are not so much different, as seen in the right panels of Fig. 2, the scattering length in the D-D * channel is larger than that in theK-D * channel. This can be attributed to the smaller kinetic energy of D in comparison toK due to a heavier charm quark. A similar tendency has also been reported in studies of phenomenological models (see e.g. [11]).
Although we find a good evidence of a sizable attraction in the I = 0 channel at m π = 410 ∼ 700 MeV, the existence of a bound or resonant

Summary
In order to clarify the possible existence of charmed tetraquark states (T cc and T cs ), we have studied the S-wave meson-meson interactions in several I = 0 and I = 1 channels (D-D,K-D, D-D * andK-D * ), using (2+1)-flavor full QCD gauge configurations generated at m π = 410 ∼ 700 MeV. For the charm quark, we have employed the relativistic heavy-quark action to take into account its dynamics on the lattice.
S-wave meson-meson potentials are extracted from Nambu-Bethe-Salpeter wave functions using the HAL QCD method. Potentials are then used to calculate scattering phase shifts and scattering lengths. S-wave meson-meson interactions in the I = 1 channels are found to be repulsive and insensitive to the pion mass in the region we explored, so that tetraquark bound states are unlikely to be formed even at the physical pion mass. On the other hand, the S-wave interactions in the I = 0 channels show attractions in thē K-D,K-D * and D-D * channels, which are qualitatively consistent with the phenomenological diquark picture. S-wave scattering phase shifts in these attractive channels indicate, however, that no bound states or resonances are formed at the pion masses used in this study, m π = 410 − 700 MeV, though attractions become more prominent as the pion mass decreases, particularly in the I = 0 D-D * channel corresponding to T cc (J P = 1 + , I = 0).
To make a definite conclusion on the fate of T cc and T cs in the real world, simulations near or at the physical point are necessary. We are planning to carry out such simulations with the PACS-CS (2+1)-flavor full QCD configurations with coupled-channel schemes [26,30].