A coupled-channel lattice study on the resonance-like structure $Z_c(3900)$

TL;DR

This lattice QCD study investigates the $Z_c(3900)$ resonance-like structure through coupled-channel scattering analysis of $D$ and $ar{D}^*$ mesons, providing insights into its nature and challenging the resonance interpretation.

Contribution

First lattice QCD analysis using coupled-channel formalism to study the $Z_c(3900)$ structure, combining finite-volume energy levels with Ross-Shaw theory to extract scattering parameters.

Findings

No peak behavior in elastic scattering cross section near threshold.

Disfavors a narrow resonance scenario for $Z_c(3900)$ beyond 3 sigma.

Scattering parameters suggest the $Z_c(3900)$ may not be a conventional resonance.

Abstract

In this exploratory study, near-threshold scattering of $D$ and $\bar{D}^*$ meson is investigated using lattice QCD with $N_f=2+1+1$ twisted mass fermion configurations. The calculation is performed within the coupled-channel L\"uscher's finite-size formalism. The study focuses on the channel with $I^G(J^{PC})=1^+(1^{+-})$ where the resonance-like structure $Z_c(3900)$ was discovered. We first identify the most relevant two channels of the problem and the lattice study is performed within the two-channel scattering model. Combined with a two-channel Ross-Shaw theory, scattering parameters are extracted from the energy levels by solving the generalized eigenvalue problem. Our results on the scattering length parameters suggest that, at the particular lattice parameters that we studied, the best fitted parameters do not correspond to a peak behavior in the elastic scattering cross section near the threshold. Furthermore, within the zero-range Ross-Shaw theory, the scenario of a narrow resonance close to the threshold is disfavored beyond $3\sigma$ level.

Tables

Table 1. Table 1: Simulation parameters in this study. The lattice spacing is about 0.0863 0.0863 0.0863 fm in physical unit while the pion mass is roughly 320 320 320 MeV.

ensemble	$\beta$	$a{\mu }_l$	$a{\mu }_{\sigma }$	$a{\mu }_{\delta }$	${\left(L/a\right)}^3\times T/a$	$N_{conf}$
A40.32	$1.9$	$0.0040$	$0.150$	$0.190$	${32}^3\times 64$	250

Table 2. Table 2: Energy eigenvalues obtained in the coupled two channel scattering.

Eigenvalues	$[t_{\min },t_{\max }]$	$E_{\alpha }$	${\chi }^2/N_{d.o.f}$
${\lambda }_0$	$[14,24]$	$1.4411\left(7\right)$	$0.6$
${\lambda }_1$	$[11,24]$	$1.5547\left(6\right)$	$1.1$
${\lambda }_2$	$[10,22]$	$1.6375\left(10\right)$	$0.7$
${\lambda }_3$	$[12,22]$	$1.6868\left(23\right)$	$0.6$
${\lambda }_4$	$[11,22]$	$1.7058\left(31\right)$	$0.3$
${\lambda }_5$	$[12,18]$	$1.7319\left(33\right)$	$0.4$
${\lambda }_6$	$[12,18]$	$1.7610\left(63\right)$	$0.4$
${\lambda }_7$	$[9,14]$	$1.8152\left(30\right)$	$0.6$

Table 3. Table 3: Parameters M 11 subscript 𝑀 11 M_{11} , M 22 subscript 𝑀 22 M_{22} , M 12 subscript 𝑀 12 M_{12} and R 𝑅 R obtained by minimizing the χ 2 superscript 𝜒 2 \chi^{2} function defined in Eq. ( 55 ). All quantities are in lattice units. The corresponding values of total χ 2 superscript 𝜒 2 \chi^{2} and the number of degrees of freedom are also listed. The errors for various parameters are obtained by doing a jackknife analysis.

No. of levels	$M_{11}$	$M_{12}$	$M_{22}$	$R$	${\chi }^2/N_{d.o.f}$
7	$-7.7\left(1.8\right)$	$3.7\left(1.1\right)$	$-3.3\left(2.8\right)$	$0.12\left(1.2\right)$	$0.57/3$
7	-7.51(0.45)	3.7(1.2)	-3.1(1.5)		0.58/4

Equations

$$\mathbf{k}=\left(\frac{2\pi }{L}\right)\mathbf{n},$$

$$E_{1\cdot 2}(\mathbf{k})=\sqrt{m_1^2+{\bar{\mathbf{k}}}^2}+\sqrt{m_2^2+{\bar{\mathbf{k}}}^2},$$

$$q^2={\bar{\mathbf{k}}}^2{L}^2/(2\pi {)}^2,$$

$$q\cot {\delta }_0(q)=\frac{1}{{\pi }^{3/2}}{\mathcal{Z}}_{00}(1;q^2),$$

$$S=\left[\begin{array}[]{cc}S_{11}&S_{12}\\ S_{12}&S_{22}\end{array}\right]=\left[\begin{array}[]{cc}\eta e^{2i\delta_{1}}&i\sqrt{1-\eta^{2}}e^{i(\delta_{1}+\delta_{2})}\\ i\sqrt{1-\eta^{2}}e^{i(\delta_{1}+\delta_{2})}&\eta e^{2i\delta_{2}}\end{array}\right]\;,$$

$$\left|\begin{array}[]{cc}\frac{{\mathcal{M}}(k^{2}_{1})+i}{{\mathcal{M}}(k^{2}_{1})-i}-S_{11}&\sqrt{\frac{k_{2}m_{2}}{k_{1}m_{1}}}S_{12}\\ \sqrt{\frac{k_{2}m_{2}}{k_{1}m_{1}}}S_{12}&\frac{{\mathcal{M}}(k^{2}_{2})+i}{{\mathcal{M}}(k^{2}_{2})-i}-S_{22}\end{array}\right|=0\;.$$

$$E=\frac{k_1^2}{2m_1}=E_T+\frac{k_2^2}{2m_2},$$

$${\sigma }_e^{(l)}$$

$${\sigma }_r^{(l)}$$

$$k\cot \delta (k)=\frac{1}{a_0}+\frac{1}{2}{r}_0{k}^2+\cdots ,$$

$$M=k^{1/2}\cdot K^{-1}\cdot k^{1/2},$$

$$k=\left(\begin{array}[]{cc}k_{1}&0\\ 0&k_{2}\end{array}\right)\;,$$

$$S=\frac{1+iK}{1-iK}.$$

$$K=\tan \delta ,$$

$$M_{ij}(E)=M_{ij}(E_0)+\frac{1}{2}{R}_i{\delta }_{ij}\left[k_i^2(E)-k_i^2(E_0)\right],$$

$$M(E)=\left(\begin{array}[]{cc}M_{11}+\frac{R_{1}}{2}[k^{2}_{1}-k^{2}_{10}]&M_{12}\\ M_{12}&M_{22}+\frac{R_{2}}{2}k^{2}_{2}\end{array}\right)\;,$$

$$K(E)=k^{1/2}\cdot \left(\sum _{\alpha =1}^n\frac{{\gamma }_{\alpha }\otimes {\gamma }_{\alpha }^T}{E-E_{\alpha }}\right)\cdot k^{1/2},$$

$$K^{-1}=T^{-1}+i,$$

$$S=1+2iT,$$

$${\sigma }_{11}=\frac{4\pi }{k_1^2}|T_{11}{|}^2.$$

$$T_{11}=\frac{1}{{\alpha }_1(E)-i},$$

$$S_{11}=e^{2i{\delta }_1^C},S_{22}=e^{2i{\delta }_2^C}.$$

$$Im({\delta }_1^C)=Im({\delta }_2^C)=\epsilon \ge 0.$$

$${\delta }_1=Re({\delta }_1^C),{\delta }_2=Re({\delta }_2^C),\eta =e^{-2\epsilon },$$

$$T_{11}=\frac{1}{\cot {\delta }_1^C-i},\cot {\delta }_1^C={\alpha }_1(E).$$

$${\alpha }_1(E)=\frac{1}{k_1(E)}\left[M_{11}(E)-\frac{M_{12}^2}{M_{22}(E)+{\kappa }_2(E)}\right],$$

$${\sigma }_e=\frac{4\pi }{k_1^2+{\left(M_{11}-\frac{M_{12}^2}{M_{22}+{\kappa }_2}\right)}^2},$$

$$M_{11}=\frac{M_{12}^2}{M_{22}+{\kappa }_2},$$

$${\kappa }_2={\kappa }_{2c}\equiv \frac{M_{12}^2}{M_{11}}-M_{22}.$$

$$\Delta \equiv M_{11}{M}_{22}-M_{12}^2.$$

$${\kappa }_{2c}=-\frac{\Delta }{M_{11}}.$$

Full text

I Introduction

In the past decade or so, various exotic hadronic resonance-like structures have been witnessed by numerous experimental groups. These structures, due to their unknown nature, are generally called $XYZ$ particles. More interesting ones are the charged structures that have been discovered in both charm and bottom sectors. These structures necessarily bear a four valence quark structure $\bar{Q}q\bar{q}^{\prime}Q$ with $Q$ being a heavy-flavor quark while $q$ and $q^{\prime}$ being the light-flavored quark. For different light flavors, these are charged objects. Another interesting feature is that they tend to appear close to the threshold of two heavy mesons with valence structure $\bar{Q}q$ and $\bar{q}^{\prime}Q$. The physical nature of these structures have been contemplated and discussed in many phenomenological studies. For example, they could be shallow bound states of the two mesons due to residual color interactions or some genuine tetraquark objects. However, even after so many phenomenological studies, the nature of many of these states remains obscure. A typical example is the structure $Z_{c}(3900)$, which will be the main topic of this paper. It was first discovered by BESIII Ablikim et al. (2013) and Belle Liu et al. (2013) and soon also verified by CLEO collaborations Xiao et al. (2013). The nature of $Z_{c}(3900)$ remains in debate. For a recent review on these matters, see e.g. Ref. Liu (2016); Guo et al. (2018). It is therefore highly desirable that non-perturbative methods like lattice QCD could provide some information on the nature of these states.

Quite contrary to many phenomenological studies, lattice studies on these states are still relatively scarce. For the state $Z_{c}(3900)$, it is readily observed that the invariant mass of the structure is close to the $D\bar{D}^{*}$ threshold, naturally suggesting a shallow molecular bound state formed by the two corresponding charmed mesons. To further investigate this possibility, the interaction between $D^{*}$ and $D$ mesons near the threshold becomes crucial.

A lattice study was performed by S. Prelovsek et al. who investigated the energy levels of the two charmed meson system in the channel where $Z_{c}$ appearing in a finite volume Prelovsek et al. (2015). They used quite a number of operators, including two-meson operators in the channel of $J/\psi\pi$, $D\bar{D}^{*}$ etc. and even tetraquark operators. However, they discovered no indication of extra new energy levels apart from the almost free scattering states of the two mesons. Taking $D\bar{D}^{*}$ as the main relevant channel, CLQCD utilized single-channel Lüscher scattering formalism Lüscher (1986a, b); Lüscher and Wolff (1990); Lüscher (1991a, b) to tackle the problem and also found slightly repulsive interaction between the two charmed mesons Chen et al. (2014, 2015). Therefore, it is also unlikely for them to form bound states. A similar study using staggered quarks also finds no clue for the existence of the state Lee et al. (2014). Thus, the above mentioned lattice studies, whether it is inspecting the energy levels alone or utilizing single-channel Lüscher approach, have found no clue for the existence of a $D\bar{D}^{*}$ bound state.

On the other hand, starting around 2015, HALQCD studied the problem using the so-called HALQCD approach Ishii et al. (2007) which is rather different from Lüscher’s adopted by other groups. They claimed that this structure can be reproduced and it is not a usual bound state or resonance, but rather, a structure formed due to strong cross-channel interactions, see Ref. Ikeda et al. (2016); Ikeda (2018) and references therein. So far, this scenario is only seen in the HALQCD approach but not within the Lüscher-type approach. In fact, the multi-channel Lüscher formula is already known Song He (2005); Liu et al. (2006); Lage et al. (2009); Döring et al. (2011, 2012). Using this multi-channel Lüscher approach, Hadron Spectrum Collaboration has successfully studied various coupled-channel scattering processes involving light mesons Dudek et al. (2014); Wilson et al. (2015a, b); Briceno et al. (2017). It is therefore tempting for us to verify/falsify the cross-channel interaction scenario suggested by HALQCD in a coupled-channel Lüscher approach. More importantly, the state $Z_{c}$ does have many coupled decay channels and if the coupled channel effects are important, then inspecting the energy levels alone or performing only a single-channel scattering study will not be adequate to understand the nature of these structures. Therefore, in this exploratory study, we aim to make a step towards the multi-channel lattice computation using the coupled-channel Lüscher approach.

It is well-known that, within single-channel Lüscher’s approach, energy levels are in one-to-one correspondence with the scattering phase. However, in a two-channel situation, the $S$-matrix is characterized by 3 parameters, all of which are functions of the scattering energy. Therefore, one needs to re-parameterize the $S$-matrix elements in terms of a number of constant parameters so as to pass from the energy levels to the scattering phases. One possible choice is to utilize the so-called $K$-matrix parameterization adopted by the Hadron Spectrum Collaboration in their studies of light meson coupled-channel scattering. In this work, however, since we are only interested in the energy region very close to the threshold, we will be using the so-called multi-channel effective range expansion developed long time ago by Ross and Shaw Ross and Shaw (1960, 1961). We will call this the Ross-Shaw theory. The difficulty with the multi-channel approach lies in the fact that, the number of parameters needed to parameterize the $S$-matrix grows quadratically fast when the number of channels is increased. Therefore, based on the experimental facts and also hints from the HALQCD study, we attempt to study this problem within the two-channel Lüscher approach. This could be viewed as the first step beyond the single-channel approximation using Lüscher’s formalism. To be more specific, we will single out two most relevant channels for $Z_{c}(3900)$: $J/\psi\pi$ and $D\bar{D}^{*}$, the first being the discovery channel for $Z_{c}(3900)$ and the second has shown to be the dominant channel by BESIII experimental data Ablikim et al. (2015). Quite similar to the single channel effective range expansion which is characterized by two real parameters, namely the scattering length $a_{0}$ and the effective range $r_{0}$, in a two-channel situation, one needs 5 real parameters to describe the so-called Ross-Shaw matrix $M$: 3 for the scattering length matrix and 2 for the effective range parameters. If one would like to go beyond two channels, these numbers go up rather quickly. For example, for the case of three channels, the scattering length matrix has 6 real parameters while the effective range will add another 3, making the total number of parameters 9 which we think is already too many to handle. Thus, in this paper, using lattice QCD within the two-channel Lüscher’s formalism combined with the Ross-Shaw effective range expansion, our aim is to check whether the cross-channel interaction scenario as suggested by the HALQCD can be realized or not.

This paper is organized as follows. In Section II, we briefly outline the computational strategies in Lüscher’s formalism and review the Ross-Shaw effective range expansion that we used to parameterize out $S$-matrix elements. In particular, we discuss the conditions that should be satisfied in order to have a narrow resonance close to the threshold. In section III, interpolating operators are introduced from which correlation matrices can be computed. We also outline how the most relevant two channels are determined from these correlation matrices. In section IV, simulation details are given and the results for the single- and two-meson systems are analyzed. By applying the two-channel Lüscher’s formula, parameters that determines the Ross-Shaw $M$-matrix are extracted and the physics implied by them are also discussed. In Section V, we will conclude with some general remarks.

II Strategies for the computation

In this section, we will outline the computational strategies for the computation described in this paper. We start by reviewing the ingredients in Lüscher’s formalism with the focus on the multi-channel version. Then we will describe the multi-channel effective range expansion developed by Ross and Shaw. Relation of the scattering cross section with respect to the parameters in Ross and Shaw theory are also outlined which will help us understand the meanings of these parameters.

II.1 Multi-channel Lüscher formula

Within the original single-channel Lüscher formalism, the exact energy eigenvalue of a two-particle system in a finite box of size $L$ is related to the elastic scattering phase of the two particles in the infinite volume. For simplicity, we will only consider the center of mass frame with periodic boundary conditions applied in all three spatial directions. Consider two interacting bosonic particles, which will be called mesons in the following, with mass $m_{1}$ and $m_{2}$ enclosed in a cubic box of size $L$. The spatial momentum ${\mathbf{k}}$ of any particle is quantized according to:

[TABLE]

with ${\bf n}$ being a three-dimensional integer. The exact energy of the two-particle system in this finite volume is denoted as: $E_{1\cdot 2}({\mathbf{k}})$. This can be obtained from lattice simulations from appropriate correlation functions. Defining a variable $\bar{{\mathbf{k}}}^{2}$ via:

[TABLE]

which would be the energy of two freely moving particles in infinite volume with mass $m_{1}$ and $m_{2}$ bearing three-momenta ${\mathbf{k}}$ and $-{\mathbf{k}}$, respectively. It is also convenient to further define a variable $q^{2}$ as:

[TABLE]

which differs from ${\bf n}^{2}$ due to the interaction between the two particles. Single-channel Lüscher’s formula gives us a direct relation of $q^{2}$ and the elastic scattering phase shift $\tan\delta(q)$ in the infinite volume: Lüscher (1991a)

[TABLE]

where ${\mathcal{Z}}_{00}(1;q^{2})$ is the zeta-function which can be evaluated numerically once its argument $q^{2}$ is given. This relation can also address the issue of bound states which is related to the phase shift $\delta(k)$ analytically continued below the threshold, see e.g. Ref. Lüscher (1991a).

For the two-channel case, the $S$-matrix now becomes a $2\times 2$ matrix in channel space. For example, for strong interaction, the $S$-matrix is usually expressed as,

[TABLE]

where $\delta_{1}$ and $\delta_{2}$ are scattering phases in channel 1 and 2 respectively and $\eta\in[0,1]$ is called the inelasticity parameter. Note that, all three parameters, $\delta_{1}$, $\delta_{2}$ and $\eta$ are functions of the energy. It is known that below the threshold, $\eta=1$ so that the coupling between the two channels are turned off kinematically.

The two-channel Lüscher formula now takes the form,

[TABLE]

The function ${\mathcal{M}}$ involves zeta-functions and the arguments $k^{2}_{1}$ and $k^{2}_{2}$ in this formula are related to the energy of the two-particle system via,

[TABLE]

with $E_{T}$ being the threshold energy. 111Although this expression is written in non-relativistic form, it is easily modified to the relativistic form. As shown by Lüscher, this is legitimate if we neglect the polarization effects (exponentially suppressed as $e^{-m_{\pi}L}$) which is what we always assume to be the case. To be specific, $E_{T}=m_{D}+m_{D^{*}}-(m_{J/\psi}+m_{\pi})$ and $m_{1}$ and $m_{2}$ should be taken as the reduced mass of the $J/\psi\pi$ and $D\bar{D}^{*}$ systems respectively.

For a given partial-wave $l$, the cross sections above the inelastic threshold consists of elastic cross section $\sigma^{(l)}_{e}$ and the reaction cross section $\sigma^{(l)}_{r}$. They are given by,

[TABLE]

where $S^{(l)}_{ij}$ being the $S$-matrix elements for partial-wave $l$ and we have utilized the unitarity condition for the $S$-matrix.

II.2 Ross-Shaw theory

As is mentioned above, since $\delta_{1}$, $\delta_{2}$ and $\eta$ are all functions of the energy, two-channel Lüscher formula (6) gives a relation among three functions. It is therefore crucial to have a parameterization of the $S$-matrix in terms of constants instead of functions and the multi-channel effective range expansion developed by Ross and Shaw Ross and Shaw (1960, 1961) serves this purpose. Here we will briefly recapitulate the major points of that theory.

In the single-channel case, this theory is just the well-know effective range expansion for low-energy elastic scattering,

[TABLE]

where $\cdots$ designates higher order terms in $k^{2}$ that vanish in the limit of $k^{2}\rightarrow 0$. Therefore, in low-energy elastic scattering, the scattering length $a_{0}$ and the effective range $r_{0}$ completely characterize the scattering process. Ross-Shaw theory simply generalize the above theory to the case of multi-channels. For that purpose, they define a matrix $M$ via

[TABLE]

where $k$ and $K$ are both matrices in channel space. The matrix $k$ is the kinematic matrix which is a diagonal matrix given by

[TABLE]

and $k_{1}$ and $k_{2}$ are related to the energy $E$ via Eq. (7). The matrix $K$ is called the $K$-matrix in scattering theory whose relation with the $S$-matrix is given by, 222$K$-matrix is hermitian so that $S$-matrix is unitary.

[TABLE]

Another useful formal expression for the matrix $K$ is

[TABLE]

where both sides are interpreted as matrices in channel space. From the above expressions, it is easily seen that $K^{-1}$ that appears in Eq. (10) is simply the matrix $\cot\delta$ and without cross-channel coupling, the $M$-matrix is also diagonal with entries $M\sim\mbox{Diag}(k_{1}\cot\delta_{1},k_{2}\cot\delta_{2})$. Thus, it is indeed a generalization of the single channel case in Eq. (9). In their original paper, Ross and Shaw showed that the $M$-matrix as function of energy $E$ can be Taylor expanded around some reference energy $E_{0}$ as,

[TABLE]

where we have explicitly written out the channel indices $i$ and $j$. The matrix $M_{ij}(E_{0})\equiv M^{(0)}_{ij}$ is a real symmetric matrix in channel space that we will call the inverse scattering length matrix; $R\equiv\mbox{Diag}(R_{1},R_{2})$ is a diagonal matrix which we shall call the effective range matrix. $k^{2}_{i}$ are the entries for the kinematic matrix defined in Eq. (11). Therefore, for two channels, there are altogether 5 parameters to describe the scattering close to some energy $E_{0}$: 3 in the inverse scattering length matrix $M^{(0)}$ and 2 in the effective range matrix $R$. As was shown in Ref. Ross and Shaw (1961) in many cases, $R_{1}\simeq R_{2}$ and we have only $4$ parameters in this description. One could further reduce the number of parameters to 3 by neglecting terms associated with effective ranges. This is called the zero-range approximation Ross and Shaw (1960). For convenience, we usually take $E_{0}$ to be the threshold of the second channel. In such a case, the two-channel Ross-Shaw $M$-matrix looks like,

[TABLE]

where $k^{2}_{10}\equiv k^{2}_{1}(E=E_{T})$.

It is understood that Ross-Shaw parameterization in Eq. (14) is equivalent to the so-called $K$-matrix parameterization with two poles. In this $K$-matrix representation, assuming there are altogether $n$ open channels, the $n\times n$ $K$-matrix is parameterized as,

[TABLE]

where $k$ is the kinematic matrix analogue of Eq. (11), the label $\alpha=1,2,\cdots,n$ designates the channels and each $\gamma_{\alpha}$ is a $1\times n$ real constant matrix (an $n$-component vector). It is shown in Ref. Ross and Shaw (1961) that this is equivalent to the effective range expansion (14) but the parameters are more flexible. In particular, $K$-matrix parameterization contains $(n^{2}+n)$ real parameters: $n^{2}$ for $n$ copies of $\gamma_{\alpha}$’s and another $n$ for the $E_{\alpha}$’s while an $n$-channel Ross-Shaw parameterization has $n(n+1)/2+n$ real parameters, $n(n-1)/2$ parameters less than the most general $K$-matrix given in Eq. (16). In this paper, we will focus on the two-channel case only.

II.3 Resonance scenario in Ross-Shaw theory

In this subsection, we investigate the possibility of a narrow peak just below the threshold of the second channel. In particular, this will be studied within the framework of two-channel Ross-Shaw theory. It turns out that this requirement will implement some constraints among the different parameters in Ross-Shaw theory. Later on, we will extract these parameters from our lattice data and therefore try to answer the question if there could exist a narrow peak in the elastic cross section.

It is convenient to inspect the resonance scenario using the so-called $T$-matrix which is continuous across the threshold. Formally, it is related to the $K$-matrix via,

[TABLE]

or equivalently as $T=K(1-iK)^{-1}$. The relation between the $S$-matrix and the $T$-matrix is given by,

[TABLE]

where both $S$ and $T$ now are $2\times 2$ matrices in channel space. Since the scattering cross section $\sigma_{ij}$ is essentially proportional to $|T_{ij}|^{2}$, in fact we have, see Eq. (8),

[TABLE]

Therefore, if we write,

[TABLE]

with $\alpha$ being real, it is seen that a resonance peak happens when $\alpha(E)=0$ and the half-width positions are at $\alpha(E)=\pm 1$ respectively. Here, of course, we have neglected the energy dependence of the kinematic factor $1/k^{2}_{1}$ in $\sigma_{11}$ which is legitimate for narrow resonances.

It is convenient to use the idea of complex phase shifts in two different channels. In this regard, one writes:

[TABLE]

In order to ensure unitarity, the imaginary part of the two complex phase shifts have to be equal and positive:

[TABLE]

The real parts of the two complex phase shifts need not have any relations. Comparing with the general representation in Eq. (5) we have

[TABLE]

and the positivity of $\varepsilon$ ensures that $\eta$ is a positive real number between zero and one.

The above equations apply when the energy is above the threshold. Below the threshold, we have practically single-channel scattering and phase shift for channel one is real and that for the second vanishes identically: $S_{11}=e^{2i\delta_{1}}$, $S_{12}=S_{21}=0$, $S_{22}=1$.

At this stage, it is important to realize a fact that the matrix $M$ which is equivalent to $\cot\delta$ could develop discontinuity at the threshold. This is understandable since usually at the new threshold, the phase shift for the newly opened channel starts from zero which is a singular point for $\cot\delta$. However, the matrix $T$, which is $e^{i\delta}\sin\delta$ is perfectly continuous across the thresholds where $\delta=0$. It is therefore more convenient to use the $T$-matrix (although parameterized by $M$-matrix) to analyze the cross sections.

Using the complex phase shift in channel one, we have

[TABLE]

We easily obtain the following formula for $\alpha_{1}(E)$,

[TABLE]

where we have substituted $-ik_{2}=\kappa_{2}$ with $\kappa_{2}>0$ for the case of a bound state in the second channel just below the threshold. 333Thus, it is readily seen that $\alpha_{1}(E)$ is indeed real.

In the zero-range approximation, meaning that we neglect the effects of the effective ranges and set both $R_{1}$ and $R_{2}$ to zero, the elastic scattering cross section reads,

[TABLE]

where $\kappa_{2}=\sqrt{2m_{2}(E_{T}-E)}$ is the binding momentum in the second channel. The function $k^{2}_{1}$ also has a mild energy dependence. The resonance occurs when the second term in the denominator exactly vanishes, giving

[TABLE]

which is equivalent to,

[TABLE]

For later convenience, we introduce the determinant of the $M$ matrix,

[TABLE]

Therefore, the value of $\kappa_{2}$ at which resonance occurs can be written as,

[TABLE]

Thus, in order to have a resonance close to the threshold the value of $\kappa_{2c}$ has to be a positive number close to zero. That means the matrix $M$ has to be singular. This also makes sense because if the matrix $M$ is singular, that means $\cot\delta$ as a matrix is singular, meaning an almost divergent scattering length, thus signaling a large scattering cross section. On the other hand, if we were to obtain a rather large value for the $M$ matrix, that means the scattering phase shift matrix itself is small, designating a small cross section.

It is also straightforward to work out the half-width points that correspond to $\alpha=\pm 1$. We call them $\kappa^{\pm}_{2}$ and they satisfy the following equation,

[TABLE]

which yields,

[TABLE]

Therefore, half-width $\Gamma$ of the peak is given by,

[TABLE]

To summarize, to characterize the narrowness of the resonance we would use the dimensionless ratio $R_{\mbox{narrow}}=\Gamma/k_{1}$ while for the closeness to the threshold, we would use the dimensionless ratio $R_{\mbox{close}}=\kappa_{2c}/k_{1}$. According to the above discussion, these two ratios read,

[TABLE]

In order to have a narrow resonance just below the threshold, both of the above ratios have to be small. It is also seen that, apart from the kinematic factor $k_{1}$, all the other quantities are determined by parameters in the $M$ matrix.

In real world, with the known information of $m_{J/\psi}=3097$MeV, $m_{\pi}=140$MeV, $m_{D}=1864$MeV and $m_{D^{*}}=2010$MeV, it is found that momentum $k_{1}=709$MeV. Therefore, for the peak of $Z_{c}(3900)$, taking values measured in Ref.: Ablikim et al. (2013), we get,

[TABLE]

both of which are small numbers. The numbers $R_{\mbox{close}}$ and $R_{\mbox{narrow}}$ may be expressed in terms of the three parameters in the $M$ matrix. To be precise, we have,

[TABLE]

Therefore, given the parameters $M_{22}$, $R_{\mbox{close}}$ and $R_{\mbox{narrow}}$ we can get the values of $M_{11}$ and $M_{12}$ which can then be compared with what we obtain from the lattice data. We will call these conditions the closeness and narrowness condition in the following.

However, since we are not having a physical pion mass in the simulation, we need to relax this constraint. For a generic pion mass, we have,

[TABLE]

In our simulation, the mass of the relevant mesons are given in the following table where all expressed in lattice units:

[TABLE]

Substitute into the expression for $k_{10}$, we get,

[TABLE]

in lattice units, corresponding to about $720$MeV in physical unit, which is close to real world value of $709$MeV. Thus, we may demand that both ratios are close to its real world values in our analysis. Note that the matrix elements in the $M$-matrix also has dimension of momentum. Therefore, it is convenient to measure everything in units of $k_{10}$. Within this unit system, the closeness and narrowness conditions read,

[TABLE]

where everything is measured in units of $k_{10}$. This is more convenient when we analyze our data later on in subsection IV.4.

III One- and two-particle operators and correlators

Single-particle and two-particle energies are measured in Monte Carlo simulations by measuring corresponding correlation functions, which are constructed from appropriate interpolating operators with definite symmetries. For the $D\bar{D}^{*}$ channel, we basically adopted the same set of operators in our previous study, see Sec. III in Ref. Chen et al. (2014). Below, we will take this channel as an example. Operators in other channels can also be constructed accordingly.

III.1 One- and two-particle operators and correlators

Let us list the interpolating operators for the charmed mesons, we utilize the following local interpolating fields in real space for $D$ mesons:

[TABLE]

together with the interpolating operator for its anti-particle ($D^{-}$): $\bar{{\mathcal{P}}}^{(d)}({\bf x},t)=[\bar{c}\gamma_{5}d]({\bf x},t)=[{\mathcal{P}}^{(d)}({\bf x},t)]^{\dagger}$. In the above equation, we have also indicated the quark flavor content of the operator in front of the definition inside the square bracket. So, for example, the operator in Eq. (41) will create a $D^{+}$ meson when acting on the QCD vacuum. Similarly, one defines ${\mathcal{P}}^{(u)}$ and $\bar{{\mathcal{P}}}^{(u)}$ with the quark fields $d({\bf x},t)$ in Eq. (41) replaced by $u({\bf x},t)$. In an analogous manner, a set of operators ${\mathcal{V}}^{(u/d)}_{i}$ are constructed for the vector charmed mesons $D^{*\pm}$ with the $\gamma_{5}$ in ${\mathcal{P}}^{(u/d)}$ replaced by $\gamma_{i}$. A single-particle state with definite three-momentum ${\mathbf{k}}$ is defined accordingly via Fourier transform, e.g.:

[TABLE]

The conjugate of the above operator is:

[TABLE]

Similar relations also hold for ${\mathcal{V}}^{(u/d)}_{i}$ and $\bar{{\mathcal{V}}}^{(u/d)}_{i}$. The interpolating operators for $J/\psi$, $\pi$, $\rho$ and $\eta_{c}$ are formed accordingly.

To form the two-particle operators, one has to consider the corresponding internal quantum numbers. Since our target state $Z^{\pm}_{c}(3900)$ state likely carries $I^{G}(J^{PC})=1^{+}(1^{+-})$, we use:

[TABLE]

Therefore, in terms of the operators defined in Eq. (41), we have used

[TABLE]

for a pair of charmed mesons with back-to-back momentum ${\mathbf{k}}$.

On the lattice, the rotational symmetry group $SO(3)$ is broken down to the octahedral group $O(Z)$. We use the following operator to create the two charmed meson state from the vacuum,

[TABLE]

where ${\mathbf{k}}_{\alpha}$ is a chosen three-momentum mode. The index $\alpha=1,\cdots,N$ with $N$ being the number of momentum modes considered in a corresponding channel. In this calculation, we have chosen $N=4$ for all channels with ${\mathbf{k}}^{2}_{\alpha}=(2\pi/L)^{2}{\bf n}^{2}_{\alpha}$, ${\bf n}^{2}_{\alpha}=0,1,2,3$. In the above equation, the notation $R\circ{\mathbf{k}}_{\alpha}$ represents the momentum obtained from ${\mathbf{k}}_{\alpha}$ by applying the operation $R$ on ${\mathbf{k}}_{\alpha}$.

In an analogous fashion, single and two-particle operators are constructed in other channels. For example, for the $J/\psi\pi$ channel, one simply replaces the operators for $D^{*}$ and $D$ by the corresponding ones for $J/\psi$ and $\pi$, respectively.

III.2 Correlation functions

One-particle correlation function, with a definite three-momentum ${\mathbf{k}}$, for the vector and pseudo-scalar charmed mesons are defined respectively as,

[TABLE]

From these correlation functions, including similar ones for other particles discussed in this study, it is straightforward to obtain the single particle energies for various lattice momenta ${\mathbf{k}}$, enabling us to check the dispersion relations for all single particles.

We now turn to more complicated two-particle correlation functions. Generally speaking, we need to evaluate a (hermitian) correlation matrix of the form:

[TABLE]

where ${\mathcal{O}}^{i}_{\alpha}(t)$ represents the two-particle operator defined in Eq. (46) in a particular channel. Two particle energies that are to be substituted into Lüscher’s formula are obtained from this correlation matrix by solving the so-called generalized eigenvalue problem (GEVP): 444We have used the matrix notation.

[TABLE]

with $\alpha=1,2,\cdots,N_{op}$ and $t>t_{0}$. The eigenvalues $\lambda_{\alpha}(t,t_{0})$ can be shown to behave like

[TABLE]

where $E_{\alpha}$ being the eigenvalue of the Hamiltonian for the system which is the quantity to be substituted into Lüscher’s formula to extract the scattering information. The parameter $t_{0}$ is tunable and one could optimize the calculation by choosing $t_{0}$ such that the correlation function is more or less dominated by the desired eigenvalues at that particular $t_{0}$ (preferring a larger $t_{0}$) with an acceptable signal to noise ratio (preferring a smaller $t_{0}$). The eigenvectors $v_{\alpha}(t,t_{0})$ are orthonormal with respect to the metric $C(t_{0})$, $v^{\dagger}_{\alpha}C(t_{0})v_{\beta}=\delta_{\alpha\beta}$ and they contain the information of the overlaps of the original operators with the eigenvectors.

III.3 sLapH smearing

To enhance the signal for the correlation matrix functions defined in the previous subsection, we have utilized the stochastic Laplacian Heavyside smearing (sLapH smearing) as discussed in Ref. Morningstar et al. (2011). Perambulators for the charm and light quarks are obtained using diluted stochastic sources. We adopted the same dilution procedure as in Ref. Helmes et al. (2015), see Sec. 2.1 in that reference for further details. The correlation functions are then constructed from these perambulators.

III.4 Singling out the most relevant two channels

There exist four relevant channels in the energy regime we are investigating, namely $J/\psi\pi$, $D\bar{D}^{*}$, $\eta_{c}\rho$ and $D^{*}\bar{D}^{*}$, with increasing thresholds. It is suggested by HALQCD collaboration that the lowest three channels had strong couplings among each other Ikeda et al. (2016); Ikeda (2018) . To verify this, we took all four channels and construct the cross-correlation matrix $\tilde{C}(t)$ as defined in Eq. (51). Within each of the four channels, in the center of mass frame, we construct the two meson operators with back-to-back three-momentum characterized by a three-dimensional integer ${\bf n}$ as in Eq. (1), starting from ${\bf n}^{2}=0$ to ${\bf n}^{2}=3$. Then, the correlation matrix is measured using the stochastically estimated perambulators obtained from the ensemble listed in Table 1. By inspecting the magnitude of the off-diagonal matrix elements relative to the diagonal ones we get a feeling about the coupling among these channels.

To visualize this, starting from the correlation matrix $C_{\alpha\beta}(t)$ defined in Eq. (48), we define a new cross-correlation matrix $\tilde{C}(t)$ at a particular temporal separation $t$ as:

[TABLE]

such that the diagonal matrix elements of $\tilde{C}(t)$ are normalized to unity. Since there are four channels and each channel has four different momenta, the matrix $\tilde{C}(t)$ is a $16\times 16$ matrix. If there were no cross-channel couplings, the matrix would be block diagonal within each channel.

The magnitude of $\tilde{C}$ is shown in Fig. 1 for $t=10$ for the case of four channels, namely $J/\psi\pi$, $D\bar{D}^{*}$, $\eta_{c}\rho$ and $D^{*}\bar{D}^{*}$. The column of the matrix is numbered from left to right according to the channels: $J/\psi\pi$, $D\bar{D}^{*}$, $\eta_{c}\rho$ and $D^{*}\bar{D}^{*}$. Within each channel, it is numbered according to increasing ${\bf n}^{2}$ for the back-to-back momenta ${\mathbf{k}}$. Similar numbering is adopted for the row of the matrix, from top to bottom. Thus the $16\times 16$ matrix is made up of $4\times 4$ block matrices, with each block of $4\times 4$ matrix representing the scattering within a particular single channel. It is seen from the figure that the coupling among channels do exist, with the $D\bar{D}^{*}$ and $J/\psi\pi$ being the most prominent ones. We remark that, this quantity $\tilde{C}$ by itself is not a physical quantity. In fact it is also dependent on the time separation $t$. Since $J/\psi\pi$ being the lightest channel among the four, its mixture with other channels will get magnified relative to other channels as $t$ increases. Nevertheless, the magnitude of the off-diagonal matrix elements of $\tilde{C}$ still offers us a qualitative description for the coupling among different channels. Since the number of parameters increases quadratically with the number of channels, in this study we are aiming at only two coupled channels. Therefore, in the following, we will focus on the two channels that are coupled most strongly, namely $D\bar{D}^{*}$ and $J/\psi\pi$. Not surprisingly, these two channels that are coupled most strongly as indicated by our simulation also coincides with what the experimental data suggested, see e.g. Ref. Ablikim et al. (2014) where the channel $D\bar{D}^{*}$ is found to be the dominant decay channel for $Z_{c}(3900)$ while the channel $J/\psi\pi$ is the discovery channel for $Z_{c}(3900)$.

IV Simulation details and results

In this paper, we have utilized $N_{f}=2+1+1$ twisted mass gauge field configurations generated by European Twisted Mass Collaboration (ETMC) at $\beta=1.9$ corresponding to a lattice spacing $a=0.0863$fm in physical units. Details of the relevant parameters are summarized in table 1.

For the measurements related to the charm sector, we have used the Osterwalder-Seiler type action with a fictitious $c^{\prime}$ quark introduced that is degenerate with $c$ Frezzotti and Rossi (2004). They form an $SU(2)$ doublet characterized by a quark mass parameter $\mu_{c}$. The up and down quark mass values are also degenerate and the parameter $\mu_{l}$ is fixed to the same value as in the sea corresponding to pion mass $m_{\pi}=320$MeV. For the charm quark, the parameter $\mu_{c}$ is fixed in such a way that the mass of $J/\psi$ on the lattice corresponds to about $2.96$GeV.

IV.1 Single meson correlators and dispersion relations

We first check the single meson dispersion relations for the mesons involved: $\pi$, $J/\psi$, $D$ and $D^{*}$. This is crucial since we need to have a well established single meson states within a finite box in order to utilize Lüscher’s formalism. We have attempted to fit the single meson correlators with various momenta using both the continuum and the lattice version of the dispersion relations:

[TABLE]

It turns out that both works fine for small enough ${\mathbf{p}}^{2}$ except that the lattice ones are better in the sense that the slope $Z_{\mbox{latt.}}$ is closer to unity than the values for $Z_{\mbox{cont.}}$. In Fig. 2, the lattice version of these dispersion relations are illustrated for the above listed mesons. Data points are obtained from the corresponding single-meson correlators while the lines representing the linear fits using Eq. (52)

IV.2 The two-particle energy levels

To extract the two-particle energy eigenvalues, we adopt the usual Lüscher-Wolff method Lüscher and Wolff (1990). For this purpose, a new matrix $\Omega(t,t_{0})$ is defined as:

[TABLE]

where $t_{0}$ is a reference time-slice. Normally one picks a $t_{0}$ such that the signal is good and stable. The energy eigenvalues for the two-particle system are then obtained by diagonalizing the matrix $\Omega(t,t_{0})$. The eigenvalues of the matrix has the usual exponential decay behavior as described by Eq. (50) and therefore the exact energy $E_{\alpha}$ can be extracted from the effective mass plateau of the eigenvalue $\lambda_{\alpha}$.

Since we are only interested in the the two-channel scenario, we focus on the correlation sub-matrix in $J\psi\pi$ and $D\bar{D}^{*}$ sector. The correlation matrix is therefore $8\times 8$ and the GEVP process yields $8$ energy levels. In order to be able to obtain a good plateau for each of these energy levels, we have tried to remove the so-called thermal-state contaminations arising from the lightest $J/\psi\pi$ state, see e.g. Ref. Helmes et al. (2015); Dudek et al. (2012). In Fig. 3 we show a typical behavior for one of the energy eigenvalue plateau plot both with and without this procedure. It is seen that, after performing the so-called thermal state removal procedure, the plateau behavior is greatly improved.

All effective mass plateaus are shown in Fig. 4 and the the final results for the $E_{\alpha}$’s and the corresponding errors, which are analyzed using jackknife method, are summarized in Table 2 together with the fitting ranges and $\chi^{2}$ per degree of freedom. It is seen that all of the energy levels develop decent plateau behavior.

IV.3 Extraction of scattering parameters

As discussed previously, using Ross-Shaw theory, we adopted a four-parameter fit for the data. To be specific, we use $M_{11}$, $M_{12}$, $M_{22}$ and $R$ as the parameters to be determined, assuming that the effective range are the same in the first and second channel. We will collectively denote these four parameters as $\{\lambda_{i}\}$, with $i=1,2,3,4$ corresponding to the above listed four parameters.

Two-channel Lüscher formula as shown in Eq. (6) can be written as,

[TABLE]

for an energy level $E$ in the finite box. The function $F$ involves various zeta-functions which can be computed numerically once the parameters are given. Thus, for a given set of parameter $\{\lambda_{i}\}$, the above formula can be viewed as an equation for $E$. In fact, one can solve for a tower of $E$ which could be compared with what is measured from lattice simulations. Therefore, following e.g. Ref. Dudek et al. (2012), we may construct a $\chi^{2}$ function as,

[TABLE]

where the summation of $\alpha$ and $\beta$ runs over all available energy levels utilized in the fit. The energy levels $E^{sol}_{\alpha}(\{\lambda_{i}\})$ are obtained by solving Eq. (54) once the parameters $\{\lambda_{i}\}$ are given. The energy levels $E^{latt}_{\alpha}$ are the corresponding ones measured from lattice data after the GEVP process. The matrix ${\mathcal{C}}$ accounts for the correlation of these quantities since they are all measured using the same set of ensemble. We estimate the covariance matrix ${\mathcal{C}}$ from our data using jackknife method.

After minimizing the function $\chi^{2}(\{\lambda_{i}\})$ numerically with respect to the four parameters, the optimal values for $\{\lambda_{i}\}$ that minimize the $\chi^{2}$ can be obtained. In our simulation, we have obtained 8 energy levels. Since the highest energy level is subject to higher states contaminations, we decided to fit for the scattering parameters both with and without this level included. It turns out that, with the highest energy level included, we keep getting rather large $\chi^{2}$ values while using the lowest $7$ levels yields a tolerable $\chi^{2}$ value. We therefore will only list the best fitted parameter values using only $7$ levels. When performing the fit, we could use the the four parameter fit with $R$ included, or in the so-called zero range approximation which sets $R=0$. It is observed that the four parameter fit yields a value of $R$ that is consistent with zero within errors. Therefore, it makes sense to perform the fit within the zero range approximation.

The covariance matrix of the fitted parameters can be estimated following a standard jackknife procedure. For the case of three parameter fit, i.e. in the zero range approximation, we obtain the inverse covariance matrix $C^{-1}$ in lattice unit as,

[TABLE]

where the column and rows of the matrix are labeled according to $(\eta_{1},\eta_{2},\eta_{3})=(M_{11},M_{12},M_{22})$. The covariance matrix $C$ itself can be easily obtained as well. If we are to transform these matrices to the unit in which $k_{10}=1$, they have to be either multiplied/divided by the value of $k^{2}_{10}$ in lattice units.

Results of the fitting procedure are tabulated in Table 3 where all dimensional quantities are in lattice units. The errors of each individual parameter is the square root of the corresponding matrix element listed the covariance matrix $C$.

IV.4 Discussion of the results

Since the presence of the effective range parameter $R$ is marginal, in the following discussion we will only focus on the case of three parameters. For later convenience, we collectively denote these parameters as $\eta=(M_{11},M_{12},M_{22})^{T}\in\mathbb{R}^{3}$ and the value of $\eta$ which minimizes the $\chi^{2}$ function will be denoted as $\eta^{*}$. We will be showing these parameters in units of $k_{10}$.

In a small neighborhood around the best fitted value for $\eta$ (in some unit), the function $\chi^{2}(\eta)$ can be parameterized as,

[TABLE]

where $w=\eta-\eta^{*}$ and the matrix $C$ is the covariance matrix for the parameters which also yields the errors (and the cross-variance) for the fitted parameters.

In any case, the matrix $C^{-1}$ is a $3\times 3$ positive-definite symmetric real matrix which can be diagonalized via some rotation matrix $R$. In fact, setting $x=R\cdot w$ and demanding that the new matrix $R^{T}C^{-1}R$ being diagonal, we have

[TABLE]

Confidence levels can then be set by using the change in the value of $\chi^{2}$ function relative to its minimum value in the parameter space. Denoting

[TABLE]

this quantity in terms of the rotated parameters $x$ is diagonal:

[TABLE]

Therefore, for a given value of $\Delta\chi^{2}(x)$, the above equation becomes a three-dimensional ellipsoid centered at the origin of $x\in\mathbb{R}^{3}$ with three half-major axis given by: $\sqrt{2\Delta\chi^{2}}\sigma_{1}$, $\sqrt{2\Delta\chi^{2}}\sigma_{2}$ and $\sqrt{2\Delta\chi^{2}}\sigma_{3}$, respectively. In terms of the original variable $w$, this is a rotated ellipsoid with the rotation characterized by the matrix $R$ which diagonalize the matrix $C^{-1}$.

It is known that, for three-parameter $\chi^{2}$ fit, the $1\sigma$, $2\sigma$ and $3\sigma$ contours 555In three dimensions, they are in fact surfaces instead of contours. But since we will be showing two-dimensional intersections, we will call them contours instead. can be set by demanding $\Delta\chi^{2}=3.53,8.02,14.2$, respectively. Thus, denoting $w=(M_{11}-M^{*}_{11},M_{12}-M^{*}_{12},M_{22}-M^{*}_{22})^{T}$, the ellipsoid,

[TABLE]

centered at $\eta^{*}=(M^{*}_{11},M^{*}_{12},M^{*}_{22})$ will enclose the $68$%, $95.4$% and $99.7$% probability as far as the parameters are concerned.

On the other hand, the closeness and narrowness condition given in Eq. (36) yields a parabola and a hyperbola on the $\eta_{1}-\eta_{2}$ plane for a given value of $\eta_{3}=M_{22}$, respectively. Therefore, if we demand that $R_{\mbox{close}}$ and $R_{\mbox{narrow}}$ to be smaller than some prescribed cut values,

[TABLE]

We can therefore check whether these conditions are supported by our $\chi^{2}$ fit.

To summarize, these conditions are, using the unit system in which $k_{10}=1$:

[TABLE]

For a given value of $\Delta\chi^{2}$, the first inequality (63) implements an ellipsoid that encloses a certain probability around the best fitted value at $\eta=\eta^{*}$; the second inequality (64) implies a region bounded by a hyperbola in the $\eta_{1}\eta_{2}$ plane, independent of the value of $\eta_{3}=M_{22}$; the third inequality, however, requires that the point $(\eta_{1},\eta_{2}$ to be in a region that is above a parabola for a given value of $\eta_{3}=M_{22}$.

Therefore, if we want to have a narrow resonance close enough to the threshold, described by two prescribed ratios: $R^{\mbox{cut}}_{\mbox{narrow}}$ and $R^{\mbox{cut}}_{\mbox{close}}$, these points have to lie within the overlapping regions of the above mentioned parabola and hyperbola. By inspecting the location of the overlapping region with respect to various constant $\Delta\chi^{2}$ contours, we are then able to set the confidence intervals for the parameters. This could be done when a particular value of $\eta_{3}=M_{22}$ is given.

In Fig. 5, the situation is illustrated at three particular values (from left to right) of $M_{22}$, namely $M_{22}=M^{*}_{22},M^{*}_{22}+\Delta M_{22},M^{*}_{22}-\Delta M_{22}$, with $\Delta M_{22}$ being the error for $M_{22}$. Here, all quantities have been converted to the unit system in which $k_{10}=1$. In each panel of the figure, the common center of the two ellipses corresponds to minimum of the $\chi^{2}$ function, i.e. the best fitted value of $(\eta^{*}_{1},\eta^{*}_{2})$ for that particular value of $\eta_{3}$. The two elliptical shaded regions around the common center correspond to the $1\sigma$ (68% probability) and $3\sigma$ (99.7% probability) contours for the parameters $(\eta_{1},\eta_{2})=(M_{11},M_{12})$. The lower left shaded region corresponds to the narrow resonance condition described by inequality (64) while the upper right shaded region corresponds to the fact that the resonance is close to the threshold, i.e. inequality (65). In this figure, the values for the two ratios have assumed their “true” physical values for $Z_{c}(3900)$, i.e. $R^{\mbox{cut}}_{\mbox{close}}=0.0211$ and $R^{\mbox{cut}}_{\mbox{narrow}}=0.065$. It is clearly seen that in this case, the overlapping regions, which are located in the top left corner in each panel, are very far away from the $3\sigma$ contours. Therefore, for this set of parameters, it is highly unlikely that the three-parameter Ross-Shaw model can describe a resonance that, like the $Z_{c}(3900)$, is both narrow and close to threshold.

Since in our simulations we do not have a physical pion mass and a physical charm quark mass, the two ratios and also the parameter $k_{10}$ could change to some values that are not their real world values. However, since the ratios are dimensionless, we do not expect them to change drastically. Similarly, as we see previously, the value of $k_{10}$ is also not very different from the physical value. Nevertheless, we still take $R^{\mbox{cut}}_{\mbox{narrow}}=R^{\mbox{cut}}_{\mbox{close}}=0.1$ and inspect the relation of the overlapping regions and the constant $\Delta\chi^{2}$ contours again. These are shown in Fig. 6. It is seen that, even at this set of parameters, the overlapping regions are still far from the $3\sigma$ contours, showing that the three-parameter Ross-Shaw model can not explain a narrow resonance close to the resonance that is described by $R^{\mbox{cut}}_{\mbox{narrow}}=R^{\mbox{cut}}_{\mbox{close}}=0.1$.

We have also tried other parameter values and the outcome is quite similar. Therefore, as far as the parameters are concerned, it seems that the three-parameter Ross-Show model cannot realize a scenario in which a narrow resonance appears close enough to the threshold. This in fact can be understood from physical arguments as follows: The reason is that, the best fitted values for the $M$ matrix elements are all very large, either in lattice unit or in $k_{10}$ unit. Recall that this matrix is related to $\cot\delta$ matrix, or the inverse scattering length matrix of the scattering, c.f. Eq. (10). Large matrix elements of $M$, if the matrix itself is nonsingular, yields large inverse scattering length, meaning a negligible scattering effect. This is the reason that the zero-range Ross-Shaw theory had a hard time to generate a narrow resonance peak near the threshold. Furthermore, if we ask why on earth that the matrix elements of $M$ turn out to be large? This is in fact implicitly hidden in our energy levels $E_{\alpha}$, see Table 2. All energy levels we utilized in the Lüscher formula are rather close to the corresponding free two-particle energy levels. This fact in turn generates large values for the $M$ matrix elements.

However, it is still premature to draw the conclusion that the three-parameter Ross-Shaw theory cannot describe a narrow resonance close to the threshold due to the following reasons: In the above arguments we have not taken into account the systematic errors. Only statistical errors are considered and they are assumed to be normally distributed. Although we used dimensionsless quantitites in our study to bypass scale setting problems, there are still several systematic effects. The result was calculated on one nonphysical ensemble. A further study is required to perform an extrapolation to the physical point. Another systematic effects is that we have only considered two channels. Of course, one could try to add more channels e.g. the $\rho\eta_{c}$, to the discussion. For that purpose, one needs to have much more energy levels since a three-channel Ross-Shaw theory, even within the zero range approximation, needs $6$ parameters. In order to nail down these, one needs more energy levels. This is also something that could be attempted in the future.

V Conclusions

Let us now outline the main conclusions from our study:

• •

In this exploratory lattice study, we utilize the coupled channel Lüscher formula together with the Ross-Shaw theory to study the near threshold scattering of $D\bar{D}^{*}$ which is relevant for the exotic state of $Z_{c}(3900)$.

• •

We single out the most strongly coupled two channels of the problem, namely $D\bar{D}^{*}$ and $J/\psi\pi$. The fact that these two particular channels show strongest coupling is supported both by our correlation matrix estimation and by the experimental facts.

• •

Our results show that the inverse scattering length parameters $M_{11}$, $M_{12}$ and $M_{22}$ are huge in magnitude, indicating that it is unlikely to generate any resonance behavior that is both narrow and close enough to the threshold.

• •

Unlike what the HALQCD collaboration finds out, our results do not support a narrow resonance-like peak close to the threshold by taking into account the most relevant two coupled channels in the problem.

• •

However, one has to keep in mind that we have not estimated the systematic uncertainties. All error estimates are purely statistical and the systematics could be due to finite lattice spacing, non-physical pion and charm quark masses, finite volume effects, etc. which needs to be clarified in future studies.

To summarize, in this paper we present an exploratory lattice study for the coupled channel scattering near the $D\bar{D}^{*}$ threshold using coupled-channel Lüscher formalism. We utilize $2+1+1$ twisted mass fermion configurations at a lattice spacing of $0.0863$fm with the pion mass value of about $320$MeV. The most two relevant channels, namely $J\psi\pi$ and $D\bar{D}^{*}$ are studied, which are singled out from four channels by a correlation matrix analysis. To extract the scattering information, we fit our lattice results using the Ross-Shaw theory, a multi-channel generalization of the conventional effective range theory. Using our lattice data, the matrix elements of $M$ matrix are obtained together with the effective range parameter, although the latter turns out to be marginal.

Our results indicate that, it is unlikely to satisfy both narrow resonance condition and the condition that is close enough to the threshold, unless the parameters happen to reside in a small corner far away from the best fitted values in our parameter space. Keep in mind that we have only considered the statistical errors, further studies with more lattice spacings and volumes, pion masses and charm quark masses, or even more channels are needed to quantify these systematic effects. We hope that this exploratory study will shed some light on the multi-channel study of charmed meson scattering which is intimately related to $Z_{c}(3900)$.

Acknowledgments

The authors would like to thank F. K. Guo, U. Meissner, A. Rusetsky, C. Urbach for helpful discussions. This work is supported in part by the Ministry of Science and Technology of China (MSTC) under 973 project ”Systematic studies on light hadron spectroscopy”, No. 2015CB856702. It is also supported in part by the DFG and the NSFC through funds provided to the Sino-Germen CRC 110 “Symmetries and the Emergence of Structure in QCD”, DFG grant no. TRR 110 and NSFC grant No. 11621131001. This work is supported in part by the National Science Foundation of China (NSFC) under the project No. 11775229, 11875169 and by the Youth Innovation Promotion Association of CAS (2015013). LL acknowledges the support from the Key Research Program of the Chinese Academy of Sciences, Grant NO. XDPB09.