The Tetraquark Candidate Zc(3900) from Dynamical Lattice QCD Simulations

TL;DR

This study uses lattice QCD simulations to analyze the structure of the Zc(3900) tetraquark candidate, concluding it is a threshold cusp rather than a conventional resonance based on coupled-channel scattering data.

Contribution

The paper provides the first lattice QCD analysis of the Zc(3900) structure, revealing it as a threshold cusp instead of a resonance, using dynamical simulations and coupled-channel scattering.

Findings

Zc(3900) is a threshold cusp, not a conventional resonance.

Interactions are dominated by off-diagonal couplings.

Lattice QCD simulations support the cusp interpretation.

Abstract

The structure of the tetraquark candidate $Z_{c}(3900)$, which was experimentally reported in $e^+ e^-$ collisions, is studied by the s-wave meson-meson coupled-channel scattering on the lattice. The s-wave interactions among the $\pi J/\psi$, $\rho \eta_{c}$ and $D \bar{D}^{*}$ channels are derived from (2+1)-flavor dynamical QCD simulations at $m_{\pi}=410$--$700$ MeV. It is found that the interactions are dominated by the off-diagonal $\pi J/\psi$-$D \bar{D}^{*}$ and $\rho \eta_{c}$-$D \bar{D}^{*}$ couplings. With the interactions obtained, the s-wave two-body amplitudes and the pole position in the $\pi J/\psi$-$\rho \eta_{c}$-$D \bar{D}^{*}$ coupled-channel scattering are calculated. The results show that the $Z_{c}(3900)$ is not a conventional resonance but a threshold cusp. A semiphenomenological analysis with the coupled-channel interaction to the experimentally observed decay mode is also presented to confirm the conclusion.

Tables

Table 1. Table 1: Parameters of the RHQ action in our simulations.

${\kappa }_Q$	$r_s$	$\nu$	$c_B$	$c_E$
0.10959947	1.1881607	1.1450511	1.9849139	1.7819512

Table 2. Table 2: Meson masses in MeV units and the number of configurations used in our simulations.

	$m_{\pi }$	$m_{\rho }$	$m_{{\eta }_c}$	$m_{J/\psi }$	$m_D$	$m_{{\overline{D}}^{\ast }}$	$N_{\mathrm{cfg}}$
Expt.	140	775	2984	3097	1870	2007
Case I	411(1)	896(8)	2988(1)	3097(1)	1903(1)	2056(3)	450
Case II	570(1)	1000(5)	3005(1)	3118(1)	1947(1)	2101(2)	400
Case III	701(1)	1097(4)	3024(1)	3143(1)	2000(1)	2159(2)	399

Table 3. Table 3: The pole positions z pole − m D − m D ¯ ∗ subscript 𝑧 pole subscript 𝑚 𝐷 subscript 𝑚 superscript ¯ 𝐷 z_{\rm pole}-m_{D}-m_{\bar{D}^{*}} in MeV unit in different sheets with three cases for the pion mass. Case I, II and III correspond to those of Table  2 in the main text. The mean value is calculated with the potential at t = 13 𝑡 13 t=13 . The first parenthesis indicates the statistical error from lattice data, and the second parenthesis indicates the systematic error from the difference between the pole position at t = 13 𝑡 13 t=13 and that at t = 15 𝑡 15 t=15 .

	[$ttb$]	[$tbb$]	[$btb$]
I	$-146\left(112\right)\left(108\right)-i38\left(148\right)\left(32\right)$	$-177\left(116\right)\left(61\right)-i175\left(30\right)\left(22\right)$	$-369\left(129\right)\left(102\right)-i207\left(61\right)\left(20\right)$
	$-93\left(55\right)\left(21\right)-i9\left(25\right)\left(7\right)$
II	$-102\left(84\right)\left(45\right)-i14\left(11\right)\left(7\right)$	$-141\left(92\right)\left(64\right)-i151\left(149\right)\left(132\right)$	$-322\left(141\right)\left(111\right)-i114\left(96\right)\left(75\right)$
	$-59\left(67\right)\left(11\right)-i3\left(12\right)\left(1\right)$
III	$-100\left(48\right)\left(29\right)-i7\left(37\right)\left(17\right)$	$-127\left(52\right)\left(43\right)-i199\left(44\right)\left(28\right)$	$-356\left(108\right)\left(28\right)-i277\left(138\right)\left(95\right)$
	$-53\left(30\right)\left(5\right)-i2\left(11\right)\left(3\right)$

Equations

$$C^{\alpha \beta }(\vec{r},t-t_0)\equiv \frac{1}{\sqrt{Z_1^{\alpha }{Z}_2^{\alpha }}}\sum _{\vec{x}}⟨0|{\varphi }_1^{\alpha }(\vec{x}+\vec{r},t){\varphi }_2^{\alpha }(\vec{x},t){\overline{\mathcal{J}}}^{\beta }(t=t_0)|0⟩,$$

$${\psi }^{\alpha }(\vec{r};W_n)=\sum _{\vec{x}}⟨0|{\varphi }_1^{\alpha }(\vec{x}+\vec{r},t){\varphi }_2^{\alpha }(\vec{x},t)|n⟩,$$

$$C^{\alpha \beta }(\vec{r},t-t_0)=\sum _n{\psi }^{\alpha }(\vec{r};W_n)A_n^{\beta }{e}^{-W_n\left(t-t_0\right)}$$

$$\left(\frac{{\nabla }^2}{2{\mu }^{\alpha }}+\frac{{\left({\vec{p}}_n^{\alpha }\right)}^2}{2{\mu }^{\alpha }}\right){\psi }^{\alpha }(\vec{r};W_n)=0\left(\left|\vec{r}\right|>R\right),$$

$$W_n=\sqrt{m_1^{\alpha }+{\left({\vec{p}}_n^{\alpha }\right)}^2}+\sqrt{m_2^{\alpha }+{\left({\vec{p}}_n^{\alpha }\right)}^2}.$$

$$K^{\alpha }(\vec{r};W_n)\equiv \left(\frac{{\nabla }^2}{2{\mu }^{\alpha }}+\frac{{\left({\vec{p}}_n^{\alpha }\right)}^2}{2{\mu }^{\alpha }}\right){\psi }^{\alpha }(\vec{r};W_n)\left(\left|\vec{r}\right|\le R\right),$$

$$K^{\alpha }(\vec{r};W_n)=\left(E_n^{\alpha }-H_0^{\alpha }\right){\psi }^{\alpha }(\vec{r};W_n)\equiv \sum _{\beta }\int d\overrightarrow{r^{\prime }}{U}^{\alpha \beta }(\vec{r},\overrightarrow{r^{\prime }}){\psi }^{\beta }(\overrightarrow{r^{\prime }};W_n).$$

$$\displaystyle\left(\begin{array}[]{cc}V_{\mathrm{LO}}^{A,A}(\vec{r})&V_{\mathrm{LO}}^{A,B}(\vec{r})\\ V_{\mathrm{LO}}^{B,A}(\vec{r})&V_{\mathrm{LO}}^{B,B}(\vec{r})\end{array}\right)$$

$$\displaystyle~{}~{}~{}~{}~{}=\left(\begin{array}[]{cc}K^{A}(\vec{r};W_{1})&K^{A}(\vec{r};W_{2})\\ K^{B}(\vec{r};W_{1})&K^{B}(\vec{r};W_{2})\end{array}\right)\cdot\left(\begin{array}[]{cc}\psi^{A}(\vec{r};W_{1})&\psi^{A}(\vec{r};W_{2})\\ \psi^{B}(\vec{r};W_{1})&\psi^{B}(\vec{r};W_{2})\end{array}\right)^{-1}.$$

$$R^{\alpha \beta }(\vec{r},t-t_0)\equiv \frac{C^{\alpha \beta }(\vec{r},t-t_0)}{e^{-(m_1^{\alpha }+m_2^{\alpha })(t-t_0)}}=\sum _n{\psi }^{\alpha }(\vec{r};W_n)e^{-\Delta W_n^{\alpha }\left(t-t_0\right)}{A}_n^{\beta },$$

$$\left(-\frac{\partial }{\partial t}-H_0^{\alpha }\right)R^{\alpha \beta }(\vec{r},t-t_0)=\sum _{\gamma }{\Delta }^{\alpha \gamma }{V}_{\mathrm{LO}}^{\alpha \gamma }(\vec{r})R^{\gamma \beta }(\vec{r},t-t_0),$$

$$R^{\alpha \beta }(\vec{r},t-t_0;A_1^{+})\equiv \frac{1}{24}\sum _{g\in \mathcal{O}}{\chi }^{(A_1^{+})}(g)R^{\alpha \beta }(g^{-1}\vec{r},t-t_0),$$

$$\displaystyle\Delta^{\alpha\beta}V^{\alpha\beta}(\vec{r})=\sum_{\gamma}\biggl{[}\left(-\frac{\partial}{\partial t}-H_{0}\right)R(\vec{r},t-t_{0};A_{1}^{+})\biggr{]}_{\alpha\gamma}\cdot\biggl{[}R(\vec{r},t-t_{0};A_{1}^{+})\biggr{]}^{-1}_{\gamma\beta}.$$

$$S_{\mathrm{RHQ}}$$

$$D(x,y)$$

$${\overline{\mathcal{J}}}_{\pi J/\psi }(t_0)$$

$${\overline{\mathcal{J}}}_{\rho {\eta }_c}(t_0)$$

$${\overline{\mathcal{J}}}_{D{\bar{D}}^{\ast }}(t_0)$$

$$V^{\alpha \beta }({\vec{p}}_{\alpha },{\vec{p}}_{\beta })=\int \frac{d\vec{r}}{(2\pi {)}^3}{e}^{-i({\vec{p}}_{\alpha }-{\vec{p}}_{\beta })\cdot \vec{r}}f(r),$$

$$t^{\alpha \beta }({\vec{p}}_{\alpha },{\vec{p}}_{\beta };E)$$

$$S^{\alpha \beta }({\vec{p}}_{\alpha },{\vec{p}}_{\beta };E)$$

$$f^{\alpha \beta }({\vec{p}}_{\alpha },{\vec{p}}_{\beta };E)$$

$$\displaystyle\left|\pi(\vec{p}_{i})\otimes\left[A(\vec{p}_{j})\otimes B(\vec{p}_{k})\right]_{\alpha}\right\rangle=\left\{\begin{array}[]{ll}\left|\pi_{1}(\vec{p}_{1})~{}\pi_{2}(\vec{p}_{2})~{}J/\psi_{3}(\vec{p}_{3});\alpha\right\rangle\\ \frac{1}{\sqrt{2}}\left(\left|\pi_{1}(\vec{p}_{1})~{}D_{2}(\vec{p}_{2})~{}\bar{D}^{*}_{3}(\vec{p}_{3});\alpha\right\rangle+\left|D_{1}(\vec{p}_{1})~{}\pi_{2}(\vec{p}_{2})~{}\bar{D}^{*}_{3}(\vec{p}_{3});\alpha\right\rangle\right)\\ \end{array}\right.,$$

$$⟨Y(4260)|T_{Y\to \pi +\alpha }(W_3)|{\vec{p}}_i,{\vec{p}}_j,{\vec{p}}_k;\alpha ⟩$$

$$=⟨Y(4260)|\left\{\Gamma +\Gamma G_0(W_3)\sum _{i=1,2}{t}_{(i)}(W_3-E_{\pi }({\vec{p}}_i),{\vec{p}}_i)\right\}|{\vec{p}}_i,{\vec{p}}_j,{\vec{p}}_k;\alpha ⟩,$$

$$⟨{\overrightarrow{p^{\prime }}}_i,{\overrightarrow{p^{\prime }}}_j,{\overrightarrow{p^{\prime }}}_k;\alpha |G_0(W_3)|{\vec{p}}_i,{\vec{p}}_j,{\vec{p}}_k;\beta ⟩$$

$${\vec{p}}_j$$

$${\vec{p}}_k$$

$$⟨{\overrightarrow{p^{\prime }}}_i,{\overrightarrow{p^{\prime }}}_j,{\overrightarrow{p^{\prime }}}_k;\alpha |t_{(i)}(W_3-E_{\pi }({\vec{p}}_i),{\vec{p}}_i)|{\vec{p}}_i,{\vec{p}}_j,{\vec{p}}_k;\beta ⟩=\delta ({\overrightarrow{p^{\prime }}}_i-{\vec{p}}_i)t_{(i)}^{\alpha \beta }({\vec{q}}_{\alpha },{\vec{q}}_{\beta },{\vec{p}}_i;W_3),$$

$$t_{(i)}^{\alpha \beta }({\vec{q}}_{\alpha },{\vec{q}}_{\beta },{\vec{p}}_i;W_3)=V_{(i)}^{\alpha \beta }({\vec{q}}_{\alpha },{\vec{q}}_{\beta })+\sum _{\gamma =\pi J/\psi ,\rho {\eta }_c,D{\bar{D}}^{\ast }}\int d{\vec{q}}_{\gamma }\frac{V_{(i)}^{\alpha \gamma }({\vec{q}}_{\alpha },{\vec{q}}_{\gamma })t_{(i)}^{\gamma \beta }({\vec{q}}_{\gamma },{\vec{q}}_{\beta },{\vec{p}}_i;W_3)}{W_3-E_{\pi }({\vec{p}}_i)-E_{jk}({\vec{p}}_i,{\vec{q}}_{\gamma })+iϵ},$$

$$V_{(i)}^{\alpha \beta }({\vec{q}}_{\alpha },{\vec{q}}_{\beta })=\sqrt{{\mathcal{N}}^{\alpha }({\vec{q}}_{\alpha })}{V}^{\alpha \beta }({\vec{q}}_{\alpha },{\vec{q}}_{\beta })\sqrt{{\mathcal{N}}^{\beta }({\vec{q}}_{\beta })},$$

$$T_{Y\to \pi +\beta }({\vec{p}}_i,{\vec{q}}_{\beta };W_3)=C_{Y\to \pi +\beta }+\sum _{\alpha =\pi J/\psi ,D{\bar{D}}^{\ast }}{f}_{\alpha }{f}_{\beta }{C}_{Y\to \pi +\alpha }$$

$$\displaystyle~{}~{}\times\biggl{(}\int d\vec{q}_{\alpha}\frac{t_{(i)}^{\alpha\beta}(\vec{q}_{\alpha},\vec{q}_{\beta},\vec{p}_{i};W_{3})}{W_{3}-E_{i}(\vec{p}_{i})-E_{j}(\vec{p}_{j})-E_{k}(\vec{p}_{k})+i\epsilon}+\int d\vec{q}_{\alpha(j)}\frac{t_{(j)}^{\alpha\beta}(\vec{q}_{\alpha(j)},\vec{q}_{\beta(j)},\vec{p}_{j}(\vec{p}_{i},\vec{q}_{\beta});W_{3})}{W_{3}-E_{i}(\vec{p}_{i})-E_{j}(\vec{p}_{j})-E_{k}(\vec{p}_{k})+i\epsilon}\biggr{)}~{},$$

$${\Gamma }_{Y\to \pi +\alpha }(W_3)$$

$$\frac{d{\Gamma }_{Y\to \pi +\alpha }(W_3)}{d{M}_{\alpha }}=(2\pi {)}^{4}8{\pi }^2\int d({\hat{p}}_i\cdot {\hat{q}}_{\alpha })\frac{E_i({\vec{p}}_i)E_j({\vec{p}}_j)E_k({\vec{p}}_k)}{W_3}{p}_i{q}_{\alpha }{\left|T_{Y\to \pi +\alpha }({\vec{p}}_i,{\vec{q}}_{\alpha };W_3)\right|}^2.$$

Full text

1 Introduction

Understanding the exotic nature of hadron structures is one of the long-standing issues in hadron physics. Exotic hadrons are different from the conventional well-established quark-antiquark states (mesons) and three-quark states (baryons). Since the discovery of the tetraquark candidate $X(3872)$ [1], many candidates of such exotic hadrons have been experimentally reported: the pentaquark states $P_{c}^{+}(4380)$ and $P_{c}^{+}(4450)$ observed by the LHCb Collaboration [2] and the tetraquark states $Z_{c}(3900)$ reported by the BESIII [3, 4], the Belle [5], and confirmed by the CLEO-c [6] Collaborations. In particular, the resonant-like $Z_{c}(3900)$ appears as a peak in both the $\pi^{\pm}J/\psi$ and $\bar{D}D^{*}$ invariant mass spectra in the $e^{+}e^{-}\to Y(4260)\to\pi^{+}\pi^{-}J/\psi$ and $\pi D\bar{D}^{*}$ reactions: the isospin quantum number is then identified as $I^{G}=1^{+}$, and the spin is favored to be $J^{PC}=1^{+-}$. If the $Z_{c}(3900)$ is a resonance, it is clearly distinct from the conventional $c\bar{c}$ states: it couples to a charmonium, yet it is charged. So, at least four quarks, $c\bar{c}u\bar{d}$ (or its isospin partners), are involved. (See the level structure and the decay scheme in Fig.1.)

There have been discussions on the internal structure of the $Z_{c}(3900)$ using phenomenological models. The $Z_{c}(3900)$ has been interpreted as a hadro-charmonium ($\pi+J/\psi$), a compact tetraquark or a $D\bar{D}^{*}$ molecule. (e.g., Refs. [7, 8]) In these cases, the $Z_{c}(3900)$ is considered to be a resonance state. On the other hand, several models developed to analyse the experimental data indicate a kinematical threshold effect (e.g., Refs. [9, 10]). However, due to the lack of information of the coupled-channel interactions among different channels (such as $\pi J/\psi$, $\rho\eta_{c}$, and $D\bar{D}^{*}$), the predictions of those models are not well under theoretical control. On the other hand, the direct lattice QCD studies with the standard method of temporal correlations show no candidate for the $Z_{c}(3900)$ eigenstate [11, 12, 13], which indicates that the $Z_{c}(3900)$ may not be an ordinary resonance state [14]. Under these circumstances, it is most desirable to execute analyses including explicitly the coupled-channel dynamics with the first-principles QCD inputs.

In this paper, we present the full account of the first lattice QCD study to determine the nature of the $Z_{c}(3900)$ on the basis of the HAL QCD method [15, 16, 17, 18]. We consider three two-body channels below $Z_{c}(3900)$ ($\pi J/\psi$, $\rho\eta_{c}$ and $D\bar{D}^{*}$) which couple with each other. In Ref. [19], we have successfully derived the coupled-channel interactions faithful to the QCD S-matrix from the equal-time Nambu-Bethe-Salpeter (NBS) wave functions on the lattice according to the coupled-channel formulation of the HAL QCD method [20, 21, 22], and then the s-wave interaction has been used to calculate ideal scattering observables such as two-body invariant mass spectra and search for the complex poles in the scattering amplitudes to unravel the nature of the $Z_{c}(3900)$. We also present invariant mass spectra of the three-body decays $Y(4260)\to\pi\pi J/\psi$ and $\pi D\bar{D}^{*}$ using the scattering amplitudes obtained in lattice QCD, and the results are then compared with experimental data. It is noted here that the HAL QCD method is theoretically identical to the Lüscher’s method [23] to obtain scattering amplitudes. Numerical tests on the consistency between the both methods has been performed in the nonresonant $\pi\pi$ channel in Refs. [24, 25] (See also Refs. [26] and [27].). We also note that the results for two-baryon systems with hyperons [28, 29, 30, 31] obtained by the coupled-channel HAL QCD method will be compared with future experimental data.

The paper is organized as follows. In Sec. 2, we present the coupled-channel formula of the HAL QCD method to calculate the scattering amplitude in the $\pi J/\psi$-$\rho\eta_{c}$-$D\bar{D}^{*}$ coupled system. Then, the numerical setup of lattice QCD simulations are summarized in Sec. 3. The results on the coupled-channel potentials, scattering amplitude and pole position are presented in Sec. 4. In Sec. 5, we discuss the three-body decay of the $Y(4260)$ and compare our results with the experimental data. Sec. 6 is devoted to summary.

2 The coupled-channel HAL QCD method

In this section, we briefly review the coupled-channel HAL QCD method [20, 21, 22, 19] by which we extract coupled-channel energy-independent potentials in hadron-hadron scatterings. We begin with a two-hadron correlation function,

[TABLE]

where each channel is specified by $\alpha=(\pi J/\psi,\rho\eta_{c},D\bar{D}^{*})$, and $\phi^{\alpha}_{i}(\vec{y},t)$ is a local Heisenberg operator at Euclidian time $t>t_{0}$ and the spatial point $\vec{y}$ for the meson $i\ (=1,2)$ with mass $m_{i}^{\alpha}$ in channel $\alpha$. $\overline{\mathcal{J}}^{\beta}(t_{0})$ denotes a two-meson operator in channel $\beta$ located at time $t=t_{0}$. The correlation function is normalized by the wave function renormalization factor given by $Z_{i}^{\alpha}$. The NBS wave function $\psi^{\alpha}(\vec{r};W_{n})$ for each scattering state specified by $n$-th energy level on the lattice is defined by

[TABLE]

and is related to Eq. (1) as

[TABLE]

with the energy eigenvalue $W_{n}$ of the $n$-th QCD eigenstate. The coefficient $A_{n}^{\beta}\equiv\langle W_{n}|\overline{\mathcal{J}}^{\beta}(t_{0})|0\rangle$ is an overlap between the eigenstate and QCD vacuum by the insertion of the source operator $\overline{\mathcal{J}}^{\beta}(t_{0})$.

Outside a hadron-hadron interaction with the range $R$, the NBS wave function in each channel $\alpha$ satisfies the free Schrödinger equation as

[TABLE]

with $\mu^{\alpha}=m^{\alpha}_{1}m^{\alpha}_{2}/(m^{\alpha}_{1}+m^{\alpha}_{2})$ being the reduced mass in channel $\alpha$. The corresponding asymptotic momentum $\vec{p}_{n}^{\alpha}$ defined in the center-of-mass (c.m.) frame is related to the energy,

[TABLE]

On the other hand, inside the interaction, the half-off-shell T-matrix is obtained by

[TABLE]

from which we define the energy-independent coupled-channel potential matrix,

[TABLE]

with $E^{\alpha}_{n}=(\vec{p}^{\alpha}_{n})^{2}/2\mu^{\alpha}$ and $H^{\alpha}_{0}=-\nabla^{2}/2\mu^{\alpha}$ in channel $\alpha$. The energy-independent coupled-channel potential $U^{\alpha\beta}(\vec{r},\vec{r^{\prime}})$ guarantees that the S-matrix is unitary and gives the correct scattering amplitude until the new inelastic threshold opens [20, 21]. The non-locality of the coupled-channel potential $U^{\alpha\beta}(\vec{r},\vec{r^{\prime}})$ is handled in terms of the velocity expansion, $U^{\alpha\beta}(\vec{r},\vec{r^{\prime}})=(V^{\alpha\beta}_{\mathrm{LO}}(\vec{r})+V^{\alpha\beta}_{\mathrm{NLO}}(\vec{r})+\cdots)\delta(\vec{r}-\vec{r^{\prime}})$, and we extract the leading order potential $V^{\alpha\beta}_{\mathrm{LO}}(\vec{r})$ in this study.

For example, $2\times 2$ coupled-channel problems (e.g., $A\leftrightarrow B$ coupled-channel with the thresholds $m_{1}^{A}+m_{2}^{A}<m_{1}^{B}+m_{2}^{B}$), the above equation Eq. (7) is solved as

[TABLE]

This implies that at least the two levels of the NBS wave functions, $\{\psi^{A}(\vec{r};W_{n}),\psi^{B}(\vec{r};W_{n})\}_{n=1,2}$ which are connected to each scattering state with a given energy $W$ in the infinite volume if $W\geq m_{1}^{B}+m_{2}^{B}$, have to be extracted from lattice QCD simulations. In our full $3\times 3$ the $\pi J/\psi$-$\rho\eta_{c}$-$D\bar{D}^{*}$ coupled system, $\{\psi^{\pi J/\psi}(\vec{r};W_{n}),~{}\psi^{\rho\eta_{c}}(\vec{r};W_{n}),~{}\psi^{D\bar{D}^{*}}(\vec{r};W_{n})\}_{n=1,2,3}$ are needed. Optimizing the source operator, $\overline{\cal J}^{\alpha}$ $(\alpha=\pi J/\psi,~{}\rho\eta_{c},~{}D\bar{D}^{*})$, in Eq. (1), one can obtain these sets of the NBS wave functions. However, as the spatial extent $L$ gets to be large, it becomes difficult to identify the energy eigen-states in the simulations. Therefore, we employ what is called the time-dependent coupled-channel HAL QCD method to calculate the coupled-channel potential, with which the energy eigen-states need not be identified.

We introduce the hadron 4-pt correlation function $R^{\alpha\beta}$ defined by

[TABLE]

with $\Delta W^{\alpha}_{n}=W_{n}-(m^{\alpha}_{1}+m^{\alpha}_{2})$. In the nonrelativistic approximation $\Delta W^{\alpha}_{n}\simeq E^{\alpha}_{n}$, the kinetic energy in Eq. (7) is replaced with the time derivative, so that $R^{\alpha\beta}$ satisfies the time-dependent Schrödinger-type equation [17, 21],

[TABLE]

with $\Delta^{\alpha\gamma}=e^{(m_{1}^{\alpha}+m_{2}^{\alpha})(t-t_{0})}/e^{(m_{1}^{\gamma}+m_{2}^{\gamma})(t-t_{0})}$. The nonrelativistic approximation can be removed, if the higher order time derivative terms associated with relativistic corrections are included, {\mathcal{O}}\bigl{(}(\partial_{t}^{2}/m^{\alpha}_{1,2})(\partial_{t}/m^{\alpha}_{1,2})^{n}\bigr{)} with $n\geq 0$, whose contributions, however, turn out to be numerically negligible in the present lattice setup with relatively large pion masses. We consider $t$ sufficiently large that the inelastic states (The lowest one is $D^{*}\bar{D}^{*}$ in the current lattice QCD setup.) become negligible in $V^{\alpha\beta}$, otherwise the inelastic channels should be taken into account explicitly in the coupled-channel scattering.

To extract the s-wave coupled-channel potential $V^{\alpha\beta}(\vec{r})=V_{\mathrm{LO}~{}(\ell=0)}^{\alpha\beta}(\vec{r})$, we project the normalized correlation function to the $A^{+}_{1}$ representation of the cubic group on the lattice,

[TABLE]

where $g\in\mathcal{O}$ are elements of the cubic group, and the associated characters of the $A_{1}^{+}$ representation are given by $\chi^{(A_{1}^{+})}(g)~{}(=1)$. Therefore, the s-wave coupled-channel is given by

[TABLE]

The systematic errors originating from higher derivative terms and the contribution from inelastic states are estimated by the $t$ dependence of scattering observables [17].

It is noted that some phenomenological parametrization to the on-shell K-matrix is necessary to approximate the energy dependence of the S-matrix given by lattice QCD simulations in the coupled-channel Lüscher’s method proposed in Refs. [32, 33, 34], while, in the HAL QCD method, the velocity expansion is employed to approximate the nonlocality of the coupled-channel potentials faithful to the S-matrix of QCD.

3 Numerical setup of LQCD simulations

We employ (2+1)-flavor QCD gauge configurations generated by the PACS-CS Collaboration [35, 36] on a $32^{3}\times 64$ lattice with the renormalization group improved gauge action at $\beta_{\rm lat}=1.90$ and the nonperturbatively $O(a)$-improved Wilson quark action at $C_{\rm SW}=1.715$. The parameters correspond to the lattice spacing $a=0.0907~{}(13)$ fm and the spatial lattice volume $L^{3}\simeq(2.9~{}{\rm fm})^{3}$. The hopping parameters are taken to be $\kappa_{ud}=0.13700,~{}0.13727,~{}0.13754$ for $u$ and $d$ quarks and $\kappa_{s}=0.13640$ for the $s$ quark. For the charm quark action, we employ the relativistic heavy quark (RHQ) action [37], which is designed to remove the leading order and next-to-leading order cutoff errors associated with the heavy charm quark mass, ${\mathcal{O}}((m_{c}a)^{n})$ and ${\mathcal{O}}((m_{c}a)^{n}(a\Lambda_{\rm QCD}))$, respectively. The RHQ action is given by

[TABLE]

In our simulation, the parameters of the RHQ action are $\kappa_{Q}$, $r_{s}$, $r_{t}$, $\nu$, $c_{B}$ and $c_{E}$, and the redundant parameter $r_{t}$ is chosen to be $1$. We take the same parameters as in Ref. [38, 39], where the relativistic dispersion relation of the 1S charmonium is reproduced. The RHQ parameters are summarized in Table 1.

The periodic boundary conditions are imposed on the three spacial directions, and the Dirichlet boundary conditions are taken for the temporal direction at $(t-t_{0})/a=\pm 32$ to avoid contaminations from the opposite propagation of mesons in time. For the sink operators, we choose the local interpolating operators, $\phi(x)=\bar{q}(x)\Gamma q(x)~{}\left(q=u,~{}d,~{}c\right)$, where $\Gamma$ denotes a $4\times 4$ matrix acting on spinor indices. We take $\Gamma=\gamma_{5}$ for pseudo-scalar mesons ($\pi$, $\eta_{c}$ and $D$) and $\Gamma=\gamma_{i}$ for vector mesons ($\rho$, $J/\psi$ and $\bar{D}^{*}$). For the source operators to create the $I=1$ two-meson states, $\pi J/\psi$, $\rho\eta_{c}$ and $D\bar{D}^{*}$, we take into account the following zero-momentum quark wall sources:

[TABLE]

For the improvement of statistics, we repeat the measurements in Eq. (1) for each gauge configuration with respect to the source position $t_{0}/a=0$ and $32$, and both forward and backward propagations are averaged to enhance the signal. Throughout this study, statistical errors are evaluated by the jackknife method.

The calculated meson masses together with their physical masses and the number of configurations $N_{\mathrm{cfg}}$ used in the simulation are summarized in Table 2. The two-meson thresholds relevant to this study are shown in Fig. 1: the $\pi\psi^{\prime}(3826)$ threshold is above the $D\bar{D}^{*}$ threshold because of the heavy pion mass in our simulation. The $\rho\to\pi\pi$ p-wave decay is not allowed with the current spatial extent $L\simeq 3$fm, so the asymptotic $\rho\eta_{c}$ is regarded as a well-defined two-body channel. Pair annihilations of charm quarks are not taken into account in the simulations.

4 Coupled-channel potential and structure of $Z_{c}(3900)$

Shown in Fig. 2 is the resulting s-wave coupled-channel potential $V^{\alpha\beta}(r)$ from Eq. (19) at time slice $(t-t_{0})/a=13$. We find that all diagonal elements of the potential, (a) $V^{\pi J/\psi,\pi J/\psi}$, (b) $V^{\rho\eta_{c},\rho\eta_{c}}$ and (c) $V^{D\bar{D}^{*},D\bar{D}^{*}}$, are weak. This observation indicates that the structure of the $Z_{c}(3900)$ is neither a simple $\pi J/\psi$ nor $D\bar{D}^{*}$ molecule. Among the off-diagonal elements of the potential, we find that the $\pi J/\psi$-$\rho\eta_{c}$ coupling in Fig. 2 (d) is also weak: this is explained by the heavy-quark spin symmetry, which tells us that the spin flip amplitudes of the charm quark are suppressed by ${\cal O}(1/m_{c})$. On the other hand, the couplings, (e) the $\pi J/\psi$-$D\bar{D}^{*}$ and (f) the $\rho\eta_{c}$-$D\bar{D}^{*}$, are both non-trivially strong: they correspond to the rearrangement of quarks between the hidden charm sector and the open charm sector. We will show later that these off-diagonal potentials play an important role to understand the structure of the $Z_{c}(3900)$.

We observe little time-slice dependence in $(t-t_{0})/a=11$–$15$ on the potential $V^{\alpha\beta}$, as shown in Fig. 3: this implies that contributions from the inelastic $D^{*}\bar{D}^{*}$ scattering states to $V^{\alpha\beta}$ are negligible, and the convergence of the derivative expansion is reliable.

As a first step to investigate the structure of the $Z_{c}(3900)$ on the basis of $V^{\alpha\beta}$ just obtained, we fit a three-range gauss function $f(r)=\sum_{n=1}^{3}g_{n}\exp[-a_{n}r^{2}]$ to the coupled-channel potential data. Using the fit function, we obtain the coupled-channel potential in the momentum representation ($p$-space) through the Fourier transformation:

[TABLE]

with $\vec{p}_{\alpha}$ being the relative momentum of the two-meson state in channel $\alpha$. Since the potential $V^{\alpha\beta}(\vec{p}_{\alpha},\vec{p}_{\beta})$ is the input to the Lippmann-Schwinger (LS) equation, let us consider the two-body T-matrix [40]:

[TABLE]

where $\vec{p}_{\alpha}$ ($\vec{q}_{\gamma}$) indicates the on-shell (off-shell) momentum of the two-meson state in channel $\alpha$ ($\gamma$). $E$ and $E_{\gamma}(\vec{q}_{\gamma})$ represent the scattering energy in the c.m. frame and the energy of the intermediate states in channel $\gamma$, respectively. The T-matrix is related to the S-matrix and the scattering amplitude $f^{\alpha\beta}(E)$ as follows:

[TABLE]

Shown in Fig. 4 are the invariant mass spectra of the two-body scattering which is the most ideal reaction process to understand the structure of the $Z_{c}(3900)$. Since the s-wave amplitude $f^{\alpha\beta}(E)$ is related to the cross section as $\sigma^{\alpha\beta}(E)=4\pi|f^{\alpha\beta}(E)|^{2}$, the two-body invariant mass spectra are given by $\mathrm{Im}f^{\alpha\alpha}(E)$ in the $\pi J/\psi$ (red circles), $\rho\eta_{c}$ (green triangles) and $D\bar{D}^{*}$ (blue squares) channels for (a) case I, (b) case II and (c) case III in Table 2. In Fig. 4, the inner errors are statistical only, while the outer ones are statistical and systematic errors added in quadrature: the systematic errors from the truncation of the derivative expansion are evaluated by the difference between $\mathrm{Im}f^{\alpha\alpha}$ at $t=13$ and that at $t=15$. The sudden enhancement in both the $\rho\eta_{c}$ and the $D\bar{D}^{*}$ spectra is caused by the opening of the s-wave thresholds. The peak structure in the $\pi J/\psi$ spectrum just above the $D\bar{D}^{*}$ threshold is induced by the $\pi J/\psi$-$D\bar{D}^{*}$ coupling. Indeed, if we switch off the off-diagonal components of $V^{\alpha\beta}$, the (red) circles turn into the (black) crosses without any peak structure at the $D\bar{D}^{*}$ threshold. This fact implies that the peak structure in the $\pi J/\psi$ spectrum [called the $Z_{c}(3900)$] is a typical “threshold cusp” [41, 40] due to the opening of the $s$-wave $D\bar{D}^{*}$ threshold.

To make sure that the $Z_{c}(3900)$ is not associated with a resonance structure but a threshold cusp, we examine the pole positions of the S-matrix on the complex energy plane according to the notation and procedure in Ref. [42] (see also Ref [19] and A). The complex energy is defined by $z=m_{1}^{\alpha}+m_{2}^{\alpha}+p_{\alpha}^{2}/2\mu^{\alpha}$, and the “top [$t$]” (“bottom [$b$]”) sheet corresponds to $0\leq\arg p_{\alpha}<\pi$ ($\pi\leq\arg p_{\alpha}<2\pi$) for the complex momentum in each channel ($\alpha=\pi J/\psi,\rho\eta_{c},D\bar{D}^{*}$). Among eight Riemann sheets for the present three-channel scattering, as shown in Fig. 5, relevant sheets to the scattering amplitudes are the [$ttt$], [$btt$], [$bbt$] and [$bbb$] sheet in the notation of Ref. [42]. If the $Z_{c}(3900)$ is a conventional resonance, the corresponding pole could be found just above the $D\bar{D}^{*}$ threshold with a negative imaginary part on the [$bbb$] sheet. Shown in Fig. 6 are the absolute value of the $\pi J/\psi$-$\pi J/\psi$ S-matrix on the (a) [$ttt$], (b) [$btt$], (c) [$bbt$] and (d) [$bbb$] sheet for case I in Table 2, and we find a pole with a large imaginary part only on [$bbb$] sheet (The pole location is schematically illustrated in Fig. 5.). The complex energy of the pole is $z_{\rm pole}-(m_{D}+m_{\bar{D}^{*}})=-167(94)(27)-i183(46)(19)$ MeV for case I, $-128(76)(33)-i157(32)(19)$ MeV for case II, and $-190(56)(42)-i44(27)(27)$ MeV for case III, with the first and second parentheses indicating the statistical and systematic errors, respectively. In all cases, the pole is located far from the $D\bar{D}^{*}$ threshold on the real axis, so that the scattering amplitude is hardly affected by the pole. With all analyses of the two-body coupled-channel scattering, we conclude that the $Z_{c}(3900)$ is not a conventional resonance but a threshold cusp.

5 Three-body decay of Y(4260)

In the previous section, it is found that the $Z_{c}(3900)$ is associated with the $D\bar{D}^{*}$ threshold cusp induced by the off-diagonal components of $V^{\alpha\beta}$. To make further connection between the result obtained from lattice QCD simulations and the experimentally observed structure in $\pi J/\psi$ and $D\bar{D}^{*}$ invariant mass spectra [3, 4, 5, 6], we analyze the three-body decays, $Y(4260)\to\pi\pi J/\psi$ and $\pi D\bar{D}^{*}$ as shown in Fig. 7, by taking into account the final state rescattering due to the coupled-channel potential extracted from lattice QCD simulations.

We define the three-body Fock space for the final states symmetric with respect to identical two pions as follows:

[TABLE]

Here $[A\otimes B]_{\alpha}=[\pi\otimes J/\psi]$, $[D\otimes\bar{D}^{*}]$ represents the interacting pair in the final state rescattering, whose flavor channels are specified by $\alpha$; the indices $i,j,k=1,2,3$ label the spectator particle with momentum $\vec{p}_{i,j,k}$ in the three-body channel $\pi+\alpha$ $(\alpha=\pi J/\psi,D\bar{D}^{*})$; all the momenta are defined in the $Y(4260)$ rest frame ($\vec{p}_{k}=-\vec{p}_{i}-\vec{p}_{j}$). In what follows, we denote the three-body final states in Eq. (32) as $|\vec{p}_{i},\vec{p}_{j},\vec{p}_{k};\alpha\rangle$ for simplicity. The three-body state is normalized as $\langle\vec{p^{\prime}}_{i},\vec{p^{\prime}}_{j},\vec{p^{\prime}}_{k};\alpha|\vec{p}_{i},\vec{p}_{j},\vec{p}_{k};\beta\rangle=\delta_{\alpha\beta}\delta(\vec{p^{\prime}}_{i}-\vec{p}_{i})\delta(\vec{p^{\prime}}_{j}-\vec{p}_{j})$.

The three-body amplitude relevant to the $Y(4260)\to\pi+\alpha$ decay is given by

[TABLE]

where $|Y(4260)\rangle$ denotes the $Y(4260)$ state; the mass of the $Y(4260)$ and the energy of the spectator pion are denoted by $W_{3}$ and $E_{\pi}(\vec{p}_{i})$, respectively. The operator $\Gamma$ represents the primary vertex of the $Y(4260)$ decay. Three-body free Green’s function $G_{0}(W_{3})$ is given by

[TABLE]

with $E_{i}(\vec{p}_{i})=\sqrt{m_{i}^{2}+\vec{p}_{i}^{2}}$. To derive the two-body t-matrix $t_{(i)}(W_{3}-E_{\pi}(\vec{p}_{i}),\vec{p}_{i})$ in the presence of a spectator pion $\pi_{i}$ ($i=1,2$) in Eq. (32), it is useful to introduce the relative momentum $\vec{q}_{\alpha}$ for an interacting $jk$-pair in its c.m. frame, then the momenta $\vec{p}_{j}$ and $\vec{p}_{k}$ are related to the momenta $\vec{p}_{i}$ and $\vec{q}_{\alpha}$ by a Lorentz boost,

[TABLE]

with the invariant mass of an interacting $jk$-pair, $M_{jk}(\vec{q}_{\alpha})=E_{j}(\vec{q}_{\alpha})+E_{k}(\vec{q}_{\alpha})$, and the energy of the pair, $E_{jk}(\vec{p}_{i},\vec{q}_{\alpha})=\sqrt{M_{jk}(\vec{q}_{\alpha})^{2}+\vec{p}_{i}^{2}}$. Using the spectator pion momentum $\vec{p}_{i}$ and relative momentum $\vec{q}_{\alpha}$ in channel $\alpha$, the two-body t-matrix is evaluated as follows:

[TABLE]

with

[TABLE]

where $V^{\alpha\beta}(\vec{q}_{\alpha},\vec{q}_{\beta})$ is given in Eq. (25), and ${\cal N}^{\alpha}(\vec{q}_{\alpha})=\mu_{i}^{\alpha}/\mu_{i}^{\alpha}(\vec{q}_{\alpha})$ is a conventional factor to account for the relativistic kinematics; $\mu_{i}^{\alpha}=m^{\alpha}_{j}m^{\alpha}_{k}/(m^{\alpha}_{j}+m^{\alpha}_{k})$ and $\mu_{i}^{\alpha}(\vec{q}_{\alpha})=E_{j}(\vec{q}_{\alpha})E_{k}(\vec{q}_{\alpha})/M_{jk}(\vec{q}_{\alpha})$ are the reduced mass and the reduced energy of particles $j$ and $k$, respectively.

Modeling the primary vertex by complex constants, $C_{Y\rightarrow\pi+\alpha}\equiv\langle Y(4260)|\Gamma|\vec{p}_{i},\vec{p}_{j},\vec{p}_{k};\alpha\rangle$ ($\alpha=\pi J/\psi,D\bar{D}^{*}$), the three-body amplitudes for $T_{Y\rightarrow\pi+\beta}$ ($\beta=\pi J/\psi,D\bar{D}^{*}$) are now given by

[TABLE]

with $(i,j)=(1,2)$ or $(2,1)$ and $(f_{\pi J/\psi},f_{D\bar{D}^{*}})=(1,1/\sqrt{2})$. The relative momentum of an interacting $ik$ pair is denoted by $\vec{q}_{\alpha(j)}$. With the three-body amplitude, the decay rate $\Gamma_{Y\to\pi+\alpha}$ is give by

[TABLE]

where $M_{\alpha}=M_{jk}(\vec{q}_{\alpha})$ is an invariant mass in channel $\alpha$. Therefore, the invariant mass spectrum is found as

[TABLE]

We note that in Eq. (41) all three particles are assumed to be in S-state. It is also noted that, in Ref. [19], we have employed the non-relativistic kinematics, but the relativistic one is used in the present study.

In the analysis of the decay rate, we employ the physical hadron masses to have the same phase space as the experiments, while $t_{(i,j)}^{\alpha\beta}$ is taken into account from the lattice data for case I. The three-body energy is fixed at the $Y(4260)$ energy ($W_{3}=4260$ MeV), and the complex couplings $C_{Y\rightarrow\pi+\alpha}$ are fitted to the $\pi J/\psi$ and $D\bar{D}^{*}$ invariant mass spectra obtained by BESIII Collaboration shown in Fig. 8 [3, 4]. For the $D\bar{D}^{*}$ invariant mass spectrum, recently updated double $D$-tag data is used [4]. Since the experimental data are in the arbitrary scale, we first focus only on the line shapes of the invariant mass spectra. In this case, we have two real parameters, $R\equiv|C_{Y\rightarrow\pi(D\bar{D}^{*})}/C_{Y\rightarrow\pi(\pi J/\psi)}|$ and $\theta\equiv\arg(C_{Y\rightarrow\pi(D\bar{D}^{*})}/C_{Y\rightarrow\pi(\pi J/\psi)})$, and the best fit values are $R=5.20(17)$ and $\theta=-11.4(3)$ degree, where we consider the statistical error only. To compare with the raw data of the experiment, we find $\mathcal{N}_{\pi J/\psi}|C_{Y\to\pi(\pi J/\psi)}|=1.59(78)\times 10^{-3}$ MeV*-2* for the $\pi J/\psi$ and $\mathcal{N}_{D\bar{D}^{*}}|C_{Y\to\pi(\pi J/\psi)}|=1.04(41)\times 10^{-3}$ MeV*-2* for the $D\bar{D}^{*}$ invariant mass distributions, where $\mathcal{N}_{\alpha}$ are the normalization factors to the raw data.

Resulting invariant mass spectra are shown in Fig. 8 where the shaded bands denote the statistical errors: we observe that the invariant mass spectra calculated with the coupled-channel potential $V^{\alpha\beta}$ well reproduce the peak structures just above the $D\bar{D}^{*}$ threshold. In both Figs. 8 (a) and 8 (b), the peak positions around $3.9$ GeV are denoted by vertical arrows. We also find that the reflection peak in Fig. 8 (a) due to the symmetrization of two identical pions. If we turn off the off-diagonal components of $V^{\alpha\beta}$ with the same constants $C_{Y\rightarrow\pi+\alpha}$, we obtain the results shown by the dashed lines, where the lines are normalized to the results obtained from the full calculations at $4$ GeV. The peak structures around $3.9$ GeV can not be reproduced without the off-diagonal components.

6 Summary

We have studied the $\pi J/\psi$-$\rho\eta_{c}$-$D\bar{D}^{*}$ coupled-channel interactions using (2+1)-flavor full QCD gauge configurations, in order to unravel the structure of the tetraquark candidate $Z_{c}(3900)$. Thanks to the HAL QCD method, we obtain the full coupled-channel potential $V^{\alpha\beta}$, whose diagonal components are all small. This indicates that the $Z_{c}(3900)$ cannot be a simple hadro-charmonium or $D\bar{D}^{*}$ molecular state.

Also, we have found the transition potential between $\pi J/\psi$ and $D\bar{D}^{*}$ is strong, which indicates that the $Z_{c}(3900)$ can be explained as a threshold cusp. To confirm this, we calculate the invariant mass spectra and pole positions associated with the coupled-channel two-body S-matrix on the basis of $V^{\alpha\beta}$. The results indeed support that the peak in the $\pi J/\psi$ invariant mass spectrum is not associated with a conventional resonance but is a threshold cusp induced by the strong $\pi J/\psi$-$D\bar{D}^{*}$ coupling. To further strengthen our conclusion, we have made a semiphenomenological analysis of the three-body decay of the $Y(4260)$, and find that the experimentally observed peak structures just above the $D\bar{D}^{*}$ threshold are well reproduced in the $Y(4260)\to\pi\pi J/\psi$ and the $Y(4260)\to\pi D\bar{D}^{*}$ decays.

To make a definite conclusion on the structure of the $Z_{c}(3900)$ in the real world, we plan to carry out full QCD simulations near the physical point. It is also an interesting future problem to study the structure of pentaquark candidates $P_{c}^{+}(4380)$ and $P_{c}^{+}(4450)$ on the basis of the coupled-channel HAL QCD method.

The author thanks all the member of the HAL QCD Collaboration for discussion. The author is also grateful to ILDG/JLDG [43] for providing us with full QCD gauge configurations used in this study, and to Doctor C.Z. Yuan for providing us with BESIII experimental data. Numerical calculations were carried out on NEC-SX9 at RCNP in Osaka University and SR16000 at YITP in Kyoto University. This study is supported in part by JSPS KAKENHI Grants Number JP17K14287, and by MEXT as “Priority Issue on Post-K computer” (Elucidation of the Fundamental Laws and Evolution of the Universe) and SPIRE (Strategic Program for Innovative REsearch).