Zc ( 3900 ) from coupled-channel scattering-how to extract hadronic interactions ? -

We present recent progress of lattice QCD studies on hadronic interactions which play a crucial role to understand the properties of atomic nuclei and hadron resonances. There are two methods, the plateau method (or the direct method) and the HAL QCD method, to study the hadronic interactions. In the plateau method, the determination of a ground state energy from the temporal correlation functions of multi-hadron systems is a key to reliably extract the physical observables. It turns out that, due to the contamination of excited elastic scattering states nearby, one can easily be misled by a fake plateau into extracting the ground state energy. We introduce a consistency check (sanity check) which can rule out obviously false results obtained from a fake plateau, and find that none of the results obtained at the moment for two-baryon systems in the plateau method pass the test. On the other hand, the HAL QCD method is free from the fake-plateau problem. We investigate the systematic uncertainties of the HAL QCD method, which are found to be well controlled. On the basis of the HAL QCD method, the structure of the tetraquark candidate Zc(3900), which was experimentally reported in e+e− collisions, is studied by the s-wave two-meson coupled-channel scattering. The results show that the Zc(3900) is not a conventional resonance but a threshold cusp. A semi-phenomenological analysis with the coupled-channel interaction to the experimentally observed decay mode is also presented to confirm the conclusion.


Introduction
One of the major goals in hadron physics is to reveal and understand various properties of atomic nuclei from quantum chromodynamics (QCD).The hadronic interactions derived from QCD become direct inputs to the ab initio many-body methods for nuclei and are expected to explain all of nuclear and hyper-nuclear systems [1,2].The hadronic interactions also play a crucial role for the studies of the equation of state [3], which is a key quantity to understand the physics of compact stars and their binary mergers [4,5].Another goal in hadron physics is to understand the properties of hadron resonances embedded into hadron scattering states.Many experiments so far have reported the candidates of the exotic hadrons, which are different from the conventional two-quark mesons and three-quark baryons, especially in the charm quark sector [6].Since the energy (mass and width) of the hadron resonances is defined as the pole of QCD S-matrix, the hadronic interactions faithful to QCD S-matrix again play an important role for a definite conclusion about the existence of the exotic candidates.For these purposes, the extraction of the hadron interactions faithful to QCD S-matrix from the first-principles lattice QCD (LQCD) calculations is an important issue.
In recent years, the two theoretically equivalent LQCD methods to extract the hadronic interactions have been employed: the Luscher's method [7] and the HAL QCD method [8,9,10,11].In both methods, the key quantity to start with is the Nambu-Bethe-Salpeter (NBS) wave function ψ(r; W n ) defined by where φ 1,2 (r, t) is a hadronic Heisenberg operator at position (r, t), and |W n represents the n-th QCD eigenstate with eigen-energy W n ; Z 1,2 is a wave function renormalization constant.The NBS wave function is shown to satisfy the Helmholtz equation (the free Schrödinger equation) outside the range of the hadronic interactions R, where W n = 2 m 2 + k 2 n denotes the relativistic energy of two identical particles with the mass m.One can show that the NBS wave function defined in Eq. ( 1) has the information on the scattering phase shifts or equivalently the S-matrix through the Helmholtz equation in Eq. (2) [7,9,11].In the Lüscher's method, the phase shifts are obtained from the energy eigenvalue W n measured in the LQCD simulations, while in the HAL QCD method, the energy-independent potential U(r, r ) defined from the NBS wave function is extracted: The HAL QCD potential U(r, r ), by definition, becomes faithful to QCD S-matrix.With this potential, the observables such as phase shifts and binding energies are calculated by solving the Schrödinger-type equation in infinite volume.Extensive applications of the HAL QCD method have been done for various multi-hadron systems at heavy quark masses [12,13,14,15,16,17,18,19,20,21,22] and near the physical quark mass [23,24,25,26,27,28].
Since both the methods utilize the properties of the same NBS wave function, they in principle give the same results and indeed agree quantitatively with each other in the case of I = 2 pion-pion scattering [29].However in the two-nucleon (NN) system the results show disagreement.It turns out, as we will see, that the determination of energies W n of two baryons is very difficult and almost impossible using so-called the "plateau/direct method", in which the energy eigenvalues are calculated from the temporal correlation of two baryons, due to the uncontrolled contamination of elastic excited states [30,31,32,33].
In this paper, we briefly review how the plateau method fails to extract the energies, and how the results obtained in the plateau method should be checked.Also we present how the (time-dependent) HAL QCD method succeeds in solving this problem [10,33].
Another advantage of the HAL QCD method is to enable us to carry out coupled-channel scatterings [34,35,36].We present our recent result of the tetraquark candidate Z c (3900) on the basis of the coupled-channel HAL QCD method [37,38].The Z c (3900) was experimentally reported as a peak in the π ± J/ψ invariant mass in the Y(4260) → ππJ/ψ decay [39,40,41].In LQCD, the finite volume spectrum of the two-meson states was obtained, and any signal of the existence of the exotic state was not found [42,43,44].To understand the nature of the Z c (3900), since the Z c (3900) may be embedded into the coupled πJ/ψ-ρη c -DD * scattering states, it is most important to extract the coupled-channel interactions and to find the pole of the S-matrix on the complex energy plane.We extract the s-wave coupled-channel meson-meson potential relevant to the Z c (3900), and scattering observables such as invariant mass spectra and the pole position of the S-matrix are calculated.Also, we calculate the decay rate of the Y(4260) → ππJ/ψ reaction and compare the result with the experiments.
The paper is organized as follows.In Sect.2, we discuss the fake-plateau problem.Then, how the HAL QCD method solves the problem is presented in Sect.3. Applying the HAL QCD method to the coupled-channel scattering, the nature of the tetraquark candidate Z c (3900) is discussed in Sect. 4. Sect. 5 is devoted to summary.

Signal-to-noise problem in multi-baryon systems
We discuss the difficulty to determine the ground state energy of multi-hadron systems in the plateau method which has been widely employed.In a single baryon case, a mass of baryons is extracted from the temporal correlation function given by where B(t) is the Heisenberg operator of baryons, and m n and a n denote the mass and the pole residue of a n-th eigenstate of one-particle states which couples to the operator B(t), respectively.The elipsis represents contributions from two or more particle scattering states.Suppose a system we are interested in has the physical spectrum shown in figure 1 (a), where only the low-lying nucleon (N with red line) exists below the inelastic Nπ scattering states (blue lines), the nucleon mass can be extracted from the temporal correlation function C B (t), where A 1 and W 1 represent the overlap factor and eigen-energy of the lowest πN scattering state, respectively.Since there exists the energy gap of The signal-to-noise ratio for the single nucleon is then evaluated as with N conf.being the number of configuration generated in simulations [45,46].This S/N shows the nucleon mass can be reliably determined from recent LQCD simulations.
On the other hand, in the two-baryon case, the temporal correlation function is given where W n corresponds to the n-th QCD eigen-energy with the presence of two-baryon interactions.Shown in figure 1  from that of single hadrons, since there exist the NN elastic excited states (green lines) in addition to the NN ground state (red line) and inelastic excited states (blue lines).The existence of the elastic excited states makes measurements of the ground state energy much more difficult, in contrast to the determination of the nucleon mass.In a finite box with a spatial extent L, the energy difference between the ground and first excited states is about m π , so that we have to take much later time data of the temporal correlation function C BB (t) than that of C B (t).The ground state energy W 1 is extracted by with t * = δE −1 .∆E = W 1 − 2m N and δE = W 2 − W 1 are the energy shift with respect to the NN threshold and the energy separation between the ground and lowest exited states, respectively.As nuclear physics suggests, we assume the energy shift ∆E is small.With the spatial extent L = 8 fm, which is employed in the recent LQCD simulations near physical point, t * is estimated as t * 8 fm.Thus, the signal-to-noise ratio for the two-nucleon system gets exponentially worse than that for the single nucleon [45,46], so that it proves that the reliable determination of the ground state energy is very difficult or even impossible in the plateau method for the systems involving two or more baryons.

Fake-plateau problem
In the previous section, we find that taking the data at late time (t > t * ) in the multi-baryon temporal correlation function is necessary to determine the ground state energies of multi-baryon systems.However, because of the worse S/N, data at early time (t = O(1) fm t * ) are often used to calculate the ground state energies.Obviously, such early time data suffer the contamination of the excited states (elastic scattering states).In fact, plateau-like structures (mirages) in the temporal correlation functions appear in the early time data caused by the uncontrolled contamination of elastic excited states as pointed out in Ref. [30].We here briefly review this fake-plateau problem (For more details in this section, see also Ref. [30]).

BB threshold
BBπ threshold δEel ΔEBB δEinel We first demonstrate how fake plateaux appear in LQCD simulations using mock-up data.We consider the normalized correlation function given by, where ∆E BB = W 1 − 2m B is the energy shift with respect to the two-baryon threshold energy, while δE el (δE inel ) is the energy difference between the ground and first excited elastic (lowest inelastic) states (see figure 2 for the corresponding spectrum).The effective energy shift becomes so that ∆E eff BB (t) − ∆E BB shows the deviation of the energy shift from its true value.In our mock-up data, we set δE el = 50 MeV and δE inel = 500 MeV, which are the typical lowest excited energies of the elastic and inelastic scattering states with L 4 fm and m π 500 MeV, respectively.We take a small value c 1 /b 1 = 0.01, since the energy difference δE inel is to be intrinsic in QCD, and one may ideally choose good operators for baryons to suppress the inelastic contaminations.On the other hand, it is much more difficult to separate the ground state from the elastic excited state by tuning the operator.It is because the energy difference between these states are not only attributed to QCD dynamics but to the box size of the LQCD simulations.So, we adopt b 2 /b 1 = 0, ±0.1 to examine the contamination of the elastic excited state.
Shown in figure 3 (a) is the time-slice dependence of ∆E eff BB (t) − ∆E BB .If the contamination of the elastic excited state is absent (black line), the effective energy shift rapidly approaches its true value ∆E BB in t ≤ 1 fm.On the other hand, if there exists the 10% contamination of the elastic excited state in R(t), we must take much later time data (t ≥ t * 8 fm) to achieve a few MeV accuracy required in nuclear physics.However, such t * is too late to obtain a good signal due to the bad S/N in practice.
Figure 3 (b) shows ∆E eff BB (t) − ∆E BB as a function of discrete time for t ≤ 2.5 fm with the typical error by random fluctuations to R(t).The plateau-like structure is found at t 1-2 fm, although the true plateau appears t 8-10 fm as shown in figure 3 (a).This plateau-like structure, with reasonable statistical errors in recent LQCD simulations, can be misleadingly identified as the true plateau.We note that the strategy to optimize the sink and/or source operator to achieve the good overlap to the ground state [47,48] does not work in the plateau method, since the mock-up data tell us that the plateau identification cannot quantify whether the operator is really optimized or not.Therefore, unless an explicit evidence is given that the obtained plateau-like structure is not a fake, the results from the plateau method should be taken with a grain of salt.
We now examine how the fake-plateau problem shows up in actual LQCD data.As an example, let us consider the ΞΞ system in which a better signal is expected than that in the NN system.The 2 + 1 flavor full QCD gauge configurations on L 3 × T = 32 3 × 48, 40 3 × 48, 48 3 × 48, 64 3 × 64 with the lattice spacing a = 0.08995 (40) fm are employed, where the spatial extents are La = 2.9, 3.6, 4.3, 5.8 fm, respectively.The pion and Ξ masses are m π = 0.51 GeV and m Ξ = 1.46 GeV, respectively.This is the same LQCD setup as in Ref. [49] and comparable to our mock-up data.To change the overlap between the ground and the elastic excited states, we employ the smeared and wall source operators in the simulation.
In figure 4 (a) and (b), the results of the ΞΞ effective energy shifts in the 1 S 0 and 3 S 1 are shown, respectively.The plateau-like structures appear at early time t = 1 fm, but the values of the effective energy shifts depend on the choice of the source operators (smeared or wall).Therefore, at least one of (or both) plateau-like structures should be fake in each figure 4 (a) and (b), and the early time data in the plateau method do not have any predictive power about the interaction energy between two hadrons.
In the literatures with the plateau method [47,48], it has been often claimed that the results from the smeared source is more reliable since it has a better overlap to the single baryon ground state.Such a claim, however, is too naive since the most severe systematics in two-baryon systems is originated from not inelastic excited states but elastic excited states.We here show evidences that at least some of the results from the smeared source suffer from uncontrolled systematic errors.
In figure 4   calculated by the linear extrapolation in 1/L 3 .∆E ∞ in the smeared source is found to be positive, and it is not physically allowed, while the resulting ∆E ∞ in the wall source is consistent with zero.Another evidence of the uncontrolled systematics can be presented by the sink operator dependence of the effective energy shift obtained from the smeared source in the 1 S 0 channel.Employing a sinksmearing function g(r) = 1 + A exp(−Br), we vary the parameters A and B. The results are shown in figure 4 (d) and indicate the strong dependence on the sink operators, while the results with the wall source are stable in varying the parameters (see Ref. [30]).
From the above results, one can conclude that the plateau-like structures at t 1 fm obtained from the smeared source are fake.Although the results obtained from the wall source do not show inconsistency with the physical conditions (infinite volume extrapolation and sink operator independence), there still remains a possibility of the fake plateau appearing in early time LQCD data.Additional evidences are necessary to confirm that the observed plateau-like structures in the wall source are the true plateaux like the black line in figure 3 (We will confirm this in the ΞΞ ( 1 S 0 ) channel in Sect.3).
EPJ Web of Conferences 175, 01023 (2018) https://doi.org/10.1051/epjconf/201817501023Lattice 2017 Combing the results given above with the general theoretical consideration using the mock-up data, we conclude that the plateau-like structures seen at t 1 fm can be a "mirage" of the temporal correlation function, and the true plateaux appear at much later time t ≥ 8 fm, so that more sophisticated method such as the variational method [50] is at least mandatory to get the reliable results of the energy shifts.

Sanity check of LQCD data
In the previous sections, we have discussed the fake-plateau (mirage) problem in the plateau method.This observation brings about a strong doubt on the existence of the NN bound state claimed by using the plateau method in Refs.[49,51,52,53,54,55,56,57]. We here introduce a method of a consistency check (sanity check) on the basis of the Lüscher's method developed in Ref. [31], which rules out certain types of false results for multi-hadron systems1 [31,32].The Lüscher's finite volume formula [7] provides the relation between the phase shifts and energy shifts on a finite box with spatial extent L as follows: for s-wave scatterings.The scattering momentum k = |k| is related to the energy shift ∆E as ∆E = 2 √ m2 + k 2 − 2m with a scattered hadron mass m.At low energies, we find the relation between k 2 and k cot δ(k) through the effective range expansion (ERE), where a, r e and P (n) represent the scattering length, the effective range and the shape parameters, respectively.
In figure 5  shows the singular parametrization (r e m π is 1-2 order of magnitude larger than its natural value r e m π ∼ O(1) with the opposite sign).When one claims a bound state in the LQCD data, it is also worth checking if a bound state pole has the physical residue.The physical pole condition [31] is given by with k 2 b corresponding to the bound state pole position 2 .The items of the sanity check are summarized as (i) consistency between ERE k 2 >0 and ERE k 2 <0 , (ii) non-singular ERE parameters, (iii) physical residue at the bound state pole.
We perform the sanity check for all existing LQCD data for the two-nucleon in the 1 S 0 and 3 S 1 channels [49,51,52,53,54,55,56,57], together with the source operator dependence.In figure 6, we show the results of the test for (a) NPL2015 [56] and (b) YIKU2012 [49], where the authors claim the bound NN at heavy pion masses.The NPL2015 data conflict with the check points (i) consistency between ERE k 2 >0 and ERE k 2 <0 and (iii) physical residue at the bound state pole, while the YIKU2012 data are against the point (ii) non-singular ERE parameters.
The results of the test for all LQCD data are summarized in table 1.As seen in table 1, all the LQCD data using the plateau method have failed the test.The results strongly indicate that the energy shifts ∆E determined by plateaux at early time suffer from uncontrolled systematic errors from the excited state contaminations.

(a) (b)
Figure 6.k cot δ(k)/m π as a function of (k/m π ) 2 for NN in 3 S 1 channel of (a) NPL2015 [56] and (b) YIKU2012 [49].The dashed lines correspond to the Lüscher's formula for each finite volume, and the black solid lines denote the bound state condition, − −(k/m π ) 2 .The red solid lines (bands) represent the best fit (statistical and systematic errors added in quadratures) of the ERE to the data at k 2 < 0, while the blue solid line (band) in figure (a) is the ERE (its error) given by NPL2015.(a) The NPL2015 data conflict with the check points (i) consistency between ERE k 2 >0 and ERE k 2 <0 and (iii) physical residue at the bound state pole.(b) The YIKU2012 data are against the point (ii) non-singular ERE parameters.The figures are taken from Ref. [31].
Table 1.A summary of sanity checks (i) consistency between ERE k 2 >0 and ERE k 2 <0 , (ii) non-singular ERE parameters and (iii) physical residue for the bound state pole, together with the source independence of ∆E.
"No" means that the source independency/sanity check has failed, while the symbol † implies there is none or only insufficient study on the corresponding item.For more details, see Ref. [31].

The time-dependent HAL QCD method
We here review briefly the time-dependent HAL QCD method developed in Ref. [10] and show how the method solves the fake-plateau problem.We also discuss the reliability of the HAL QCD method through several tests.

Formalism
In LQCD, we start with the normalized hadron-hadron correlation function R(r, t), with A n = W n |J(t = 0)|0 .J(t = 0) stands for a source operator at t = 0 which creates two-hadron states, and φ 1,2 is an interpolating sink operator for hadrons with the corresponding wave function renormalization constant Z 1,2 .W n = 2 m 2 + k 2 n is the relativistic energy of the n-th eigenstate |W n , and E th represents an inelastic threshold energy.ψ(r, W n ) represents the NBS wave function for the elastic scattering states defined in Eq. ( 1), where the energy is restricted in W n < E th .
Consider t sufficiently large that the contributions from elastic scattering states and possible bound states remain, while those from inelastic states are negligible, thus t = O(1) fm.Then, the HAL QCD potential U is defined [9] as a solution of for all elastic eigenstates |W n .with H 0 = −∇ 2 /m and E n = k 2 n /m.The non-local but energyindependent potential U(r, r ) can be shown to exist, by explicitly constructing it as where ψ * (r; W n ) is the dual basis associated with ψ(r; W n ), and the summation over n is restricted to the elastic channel (For more details, see [9,11]).Note that while the U(r, r ) depends on the definition of the sink operator (the "scheme" in the HAL QCD method), physical observables such as phase shifts and binding energies are uniquely obtained from U(r, r ).In principle, using Eq. ( 16), the potential is extracted from R(r, t) at late time t = O(10) fm.In practice, however, R(r, t) is much noisier at late time, as demonstrated in Sect.2, so that an accurate determination of the potential in this way becomes difficult.In order to overcome this practical difficulty, the time-dependent HAL QCD method developed in [10] is employed.Since all the elastic states ψ(r, W n < E th ) share the same energy-independent potential U(r, r ) defined in Eq. ( 16), we can extract the energy-independent potential from the time-dependent Schrödinger-type equation for R(r, t) given by 1 4m The advantage of the time-dependent HAL QCD method is that one can extract the potential dictating all the elastic scattering states below the inelastic threshold from the LQCD data at t 1 fm.This is The non-locality of the potential U can be handled in terms of the velocity expansion.At the leading order of the expansion, we obtain the effective potential U(r, r ) V eff (r)δ(r − r ) as Including the next-to-leading (NLO) order term, we extract both the leading order V LO and the NLO V NLO potentials from for several R(r, t) with different source operators.

Reliability check of the HAL QCD method
We examine the reliability of the time-dependent HAL QCD method and show that systematic uncertainties, such as the source operator dependence and convergence of the velocity expansion, are under control.We also show the consistency between the HAL QCD and Lüscher's methods.On the basis of the HAL QCD method, we confirm that the mirage plateau introduced in Sect.2.2 indeed manifests as the contamination of elastic excited states (For more details, see also Ref. [33]).We employ the 2 + 1 flavor full QCD gauge configurations on L 3 × T = 40 3 × 48, 48 3 × 48, 64 3 × 64 with the lattice spacing a = 0.08995(40) fm, which are the same LQCD setup used in Sect.2.2 for the study of the ΞΞ system; The pion and Ξ masses are m π = 0.51 GeV and m Ξ = 1.46 GeV, respectively.For the source operators, we employ the smeared and wall source operators in the simulation.The choice of the sink operator corresponds to the choice of the "scheme" in the HAL QCD method.We employ the point sink operator (point-sink scheme), the standard scheme in the HAL QCD method, with the baryon operator, abc (q T a Cγ 5 q b )q c with subscripts a, b, c and q being the color indices and the quark operator, respectively.As already shown in figure 4, the energy shift ∆E from the temporal correlation function strongly depends on the choice of the source operators.Keeping this in mind, shown in figure 7 (a) is the source operator dependence of the normalized correlation function R(r, t) in Eq. ( 15) for time-slice t = 10-16.R(r, t) shows the strong source dependence, as it is expected from figure 4. R(r, t) with the smeared source (R smear (r, t)) shows strong t-dependence, while R(r, t) with the wall source (R wall (r, t)) is stable in t.This fact indicates that there exist the sizable contamination of the elastic excited states in R smeared .It is also implied that the wall source has a good overlap to the ground state, and the contamination of the elastic excited states in R wall (r, t) seems already suppressed at t 1 fm.We will quantitatively confirm this later.
Figure 7 (b) shows the effective LO potential V eff calculated from R wall (r, t) in Eq. ( 19) and the LO potential V LO obtained from both R smear (r, t) and R wall (r, t) in Eq. ( 20).The effective LO potential from the wall source V wall eff is almost identical to the LO potential V LO .

(a) (b)
Figure 8.(a) The ΞΞ( 1 S 0 ) scattering phase shift calculated from V wall eff (green), V LO (red) and V LO + V NLO ∇ 2 (blue) potentials.(b) k cot δ(k)/m π as a function of (k/m π ) 2 .The dashed lines correspond to the Lüscher's formula.Open circles are obtained from the finite volume Hamiltonian with V wall eff , while the red band shows k cot δ(k)/m π obtained from V wall eff by solving the Schrödinger equation in the infinite volume.The figures are taken from Ref. [33].
In order to check the convergence of the velocity expansion, we calculate the ΞΞ( 1 S 0 ) scattering phase shift using V wall eff , V LO and V LO + V NLO ∇ 2 potentials.In figure 8 (a), the resulting phase shifts are shown, and it is found that the ΞΞ( 1 S 0 ) system is attractive but unbound.These three potentials give the consistent results within the statistical errors at low energies, and the NLO correction appears only at high energies.The results show that systematic uncertainties coming from the velocity expansion is well under control in the point-sink scheme, and the effective LO potential and resulting scattering phase shift from the wall source are reliable at low energies in the this system.The good convergence of the velocity expansion in the point-sink scheme has been also confirmed for the NN system [12] and I = 2 ππ system [22,29] (See also Ref. [22] for the study of the smeared-sink scheme).
Next, we show the consistency between the HAL QCD and Lüscher's methods.Using the effective LO potential V wall eff , we solve the Schrödinger equation in the finite volume and obtain the energy eigenvalues ∆E(L) = 2 √ m 2 + k 2 − 2m.Through the Lüscher's finite volume formula in Eq. ( 12), we calculate k cot δ(k) at each L. The results are shown in figure 8  We finally reveal the origin of the fake plateau emerging in the plateau method in figure 4. Solving the Schrödinger equation in the finite volume with V wall eff , we get n-th eigen-functions Ψ n (r) and corresponding eigen-energies ∆E n (The reliability of these eigen-functions and eigen-energies is established in Ref. [33]).With them, the normalized correlation function R(r, t) in Eq. ( 15) is decomposed as r R smear/wall (r, t) where a smear/wall n are determined by the orthogonality of Ψ n (r).Using three low-lying eigen-functions, we reconstruct the effective energy shift given by and the results are shown in figure 9 (a).The reconstructed effective energy shift well reproduces (fake) plateau structure at t 15a, and show the true ground state appears around t ∼ 10 fm for the smeared source, while the result with the wall source is dominated by the ground state already at t ∼ 1 fm.This is understood from figure 9 15).The wall source results (red circles) show that the contamination of the first elastic excited state is less than 1%, and therefore we observe the weak time-slice dependence and the ground state dominance in R wall (r, t) at t 1 fm.On the other hand, the smeared source results show about 10% contamination of the first elastic excited state (with negative sign), and the fake plateau appears at early time t 1 fm in R smear (r, t) (see also figure 3).Note that the smallness of elastic excited state contaminations in the wall source is not guaranteed in general, so it is always essential to employ a theoretical framework which does not require the ground state saturation.
In summary, we have established the reliability of the HAL QCD method: (i) taming the elastic excited states, (ii) the convergence in the velocity expansion of non-local interaction kernels and (iii) consistency with the Lüscher's method.It is also emphasized that, by the HAL QCD method, the origin of the fake-plateau problem in the plateau method is quantitatively revealed as the uncontrolled contamination of the elastic excited states.

Nature of the Z c (3900) from coupled-channel scattering
Recently, the tetraquark candidate Z c (3900) has been reported by the BESIII [39], the Belle [40] and the CLEO-c [41] Collaborations, where a resonant-like peak appears in both the π ± J/ψ and DD * invariant mass spectra in the e + e − → Y(4260) → π + π − J/ψ and π DD * 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 obviously distinct from the conventional cc states: it couples to a charmonium, yet it is charged.So, at least four quarks, ccu d (or its isospin partners), are involved.(See the level structure and the decay scheme in figure 10.) Despite the experimental observations, the direct lattice QCD studies with the variational method of temporal correlations show no candidate for the Z c (3900) eigenstate [42,43,44], which indicates that the Z c (3900) may not be an ordinary resonance state [58].Under these circumstances, it is most desirable to execute analyses including explicitly the coupled-channel dynamics with the firstprinciples QCD inputs.
Here, we present the lattice QCD study to determine the nature of the Z c (3900) on the basis of the HAL QCD method [8,9,10,11], which is well-established as shown in Sect.3. We consider three two-body channels below Z c (3900) (πJ/ψ, ρη c and DD * ) which couple with each other.In Refs.[37,38], we have successfully derived the coupled-channel interaction faithful to the QCD Smatrix from the NBS wave function on the lattice according to the coupled-channel formulation of the HAL QCD method [34,35,36], 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) → ππJ/ψ using the scattering amplitudes obtained in lattice QCD, and the results are then compared with experimental data.

The coupled-channel HAL QCD method
We present the coupled-channel HAL QCD method [34,35,36,37] by which we extract coupledchannel energy-independent potentials in hadron-hadron scatterings.The derivation is in parallel to that in single channel channel scatterings given in Sect.3.1.The NBS wave function of channel α, ψ α (r; W n ), for each scattering state specified by n-th energy level on the lattice is defined by where each channel is specified by α = (πJ/ψ, ρη c , DD * ).Outside the hadron-hadron interaction with the range R, each NBS wave function satisfies the free Schrödinger equation as with is the reduced mass in channel α.Using Eq. ( 25), one can show that the above coupled-channel NBS wave functions have the information on the coupled-channel S-matrix [34].Inside the interaction, we define the coupled-channel potential as, Such potential is energy-independent, as explicitly shown by where ψ * β (r; W n ) is the dual basis associated with ψ β (r; W n ) in channel β and the summation over n is restricted to the energy region below the inelastic threshold energy that we do not consider in the simulations (For more details on construction of the coupled-channel potential, see Ref. [35]).The non-locality of the coupled-channel potential U αβ (r, r ) is handled in terms of the velocity expansion, , and we extract the leading order potential V αβ LO (r) in this study.
Introducing the normalized correlation function, (t = 0) denoting a two-meson operator in channel β located at time t = 0, we end up with the time-dependent Schrödinger-type equation to extract the LO coupled-channel potential [10,35]: with In the above Eq.( 29), we neglect the relativistic correction associated with the higher order time derivatives ) n with n ≥ 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 * D * in the current lattice QCD setup) become negligible in V αβ , otherwise the inelastic channels should be taken into account explicitly in the coupled-channel scattering.
EPJ Web of Conferences 175, 01023 (2018) https://doi.org/10.1051/epjconf/201817501023Lattice 2017 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 coupledchannel Lüscher's method proposed in Refs.[59,60,61], 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.

Numerical setup
In order to extract V αβ (r) from lattice QCD simulation, we employ the (2+1)-flavor QCD gauge configurations generated by the PACS-CS collaboration [62,63] on a 32 3 × 64 lattice with the renormalization group improved gauge action at β lat = 1.90 and the non-perturbatively O(a)-improved Wilson quark action at C SW = 1.715.The parameters correspond to the lattice spacing a = 0.0907( 13) fm (the spatial extent L 2.9 fm).The hopping parameters are taken to be κ ud = 0.13700, 0.13727, 0.13754 for u and d quarks and κ s = 0.13640 for the s quark.For the charm quark, we employ the relativistic heavy quark action [64] to remove the leading and next-to-leading order cutoff errors, O((m c a) n ) and O((m c a) n (aΛ QCD )), respectively [65].The calculated meson masses are listed in Table 2 together with their physical values.The two-meson thresholds relevant to our analysis are shown in figure 10: due to the heavy pion mass in our simulation, the πψ(2S ) threshold is above the DD * threshold.Also, ρ → ππ decay is not allowed with L 3fm, so that ρη c is a well defined two-body channel.The pair annihilations of charm quarks are not taken into account in the present simulations.

The coupled-channel potential and structure of the Z c (3900)
In figure 11, we show the results of the s-wave πJ/ψ-ρη c -DD * coupled-channel potentials at timeslice t = 13, where the time-slice dependence in t = 11-15 on the potentials V αβ is found to be weak: this implies that contributions from the inelastic D * D * scattering states to V αβ are negligible, and the convergence of the derivative expansion is reliable.We find that all diagonal potentials, (a) V DD * , DD * , (c) V ρη c ,ρη c and (f) V πJ/ψ,πJ/ψ , are very weak.This observation indicates that the Z c (3900) is neither simple πJ/ψ nor DD * molecule.Among the off-diagonal potentials, we find that the πJ/ψ-ρη c coupling in figure .11 (e) is also weak: this is consistent with the heavy-quark spin symmetry, which tells us that the spin flip amplitudes of the charm quark are suppressed by O(1/m c ).On the other hand, (b) the ρη c -DD * coupling and (d) the πJ/ψ-DD * coupling are both strong: they correspond to the rearrangement of quarks between the hidden charm sector and the open charm sector.
To investigate the structure of the Z c (3900) on the basis of V αβ just obtained, let us consider the two-body t-matrix, where k α (q γ ) denotes the on-shell (off-shell) momentum of the two-meson state in channel α (γ).W c.m. and E γ (q γ ) represent the scattering energy in the center-of-mass frame and the energy of the intermediate states in channel γ, respectively.Shown in figure 12 (a) are the invariant mass spectra Im f αα (W c.m. ) = −πµ α Imt αα (W c.m. ) in the πJ/ψ (red circles), ρη c (green triangles) and DD * (blue squares) channels obtained from LQCD for case I in table 2. The amplitude f αβ (W c.m. ) is related to the differential cross section as dσ αβ /dΩ = | f αβ (W c.m. )| 2 .In figure 12 (a), the shorter errors are statistical only, while the longer 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 Im f αα at t = 13 and that at t = 15.The peak structures in the ρη c and DD * spectra are caused by the opening of the s-wave thresholds.We observe the sudden enhancement in the πJ/ψ spectrum just above the DD * threshold induced by the strong πJ/ψ-DD * coupling.Indeed, if we switch off the off-diagonal components of V αβ , the red circles turn into the black crosses without showing any peak structure.This implies that the peak structure in the πJ/ψ spectrum (called Z c (3900)) is a typical "threshold cusp" [66,67] due to the opening of the s-wave DD * threshold.
Further evidence that the Z c (3900) is not associated with the resonance structure is found by the pole positions of the S-matrix on the complex energy plane.According to the notation and procedure in Ref. [68], the complex energy is defined as z = m α 1 + m α 2 + p 2 α /2µ α , and the "top [t]" ("bottom [b]") sheet corresponds to 0 ≤ arg p α < π (π ≤ arg p α < 2π) for the complex momentum p α in each channel (α = πJ/ψ, ρη c , DD * ).Among eight Riemann sheets for the present three-channel system, the most relevant one to the Z c (3900) energy region is the [bbb] sheet in the notation of [68].We find a pole  2. In figure 12 (a), the shorter error is statistical, while the longer one is statistical and systematic combined in quadrature.The figure is taken from Ref. [37].
with a large imaginary part on this sheet: z − (m D + m 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 parenthesis indicating the statistical and systematic errors, respectively.Shown in figure 12 (b) is the complex pole on the [bbb] sheet for case I.It is located far from the DD * threshold on the real axis, so that the amplitude is hardly affected by the pole.Also, on the other relevant Riemann sheets, we do not find any pole of the S-matrix relevant to the peak structure of the Z c (3900).In order to make further connection between the above results and the experimentally observed structure in πJ/ψ invariant mass spectra [39,40,41], let us now consider a semi-phenomenological analysis of the three-body decay Y(4260) → ππJ/ψ given in figure 13, by taking into account the final state rescattering due to V αβ extracted from LQCD simulations.Modeling the primary vertex by complex constants C Y→π+α (α = (πJ/ψ, DD * )), the three-body T-matrix T Y→ππJ/ψ is given by

Three-body decay of
where W 3 , E π i (k π i ) and E α (k π i , q α ) (i = 1, 2) represent the energies of the Y(4260), the spectator pion (π i ) with the momentum k π i and the interacting pairs with the relative momentum q α in channel α, respectively.It is noted that the two-body t-matrix t αβ (i) depends not only on the two-body relative momentum p π i ψ in the final state, but explicitly on the spectator pion momentum k π i (see Ref. [38]).The decay rate in the rest frame of Y( 4260) is then obtained as In the analysis, we employ the physical hadron masses in order to have the same phase space as the experiments, while t αβ is taken from the lattice data for case I.The complex couplings C Y→π+α are fitted to the BESIII data [39].The resulting decay spectrum is shown in Fig. 14, where the shaded band denotes the statistical error: we find that the coupled-channel potential V αβ well reproduces the peak structure just above the DD * threshold at 3.9 GeV.If we turn off the off-diagonal components of V αβ with the same constants C Y→π+α , we obtain the result shown by the gray dashed line, where the line is normalized to the results obtained from the full calculation at 4 GeV.The peak structure at 3.9 GeV disappears in this case.This fact shows that the key to understand the nature of the Z c (3900) is the off-diagonal coupling between the πJ/ψ and DD * channels, and the peak observed in the experiments can be reproduced by the two-body amplitude t αβ in which not a resonance pole but a threshold cusp is developed.

Summary
We have addressed the fake-plateau (mirage) problem for multi-baryon systems in the commonly used plateau method in which one tries to determine the ground state energy by the temporal correlation functions.There appear the plateau-like structures in the effective energy shifts ∆E(t) at t = O(1) fm despite the existence of the contamination of the elastic excited state, and one may be misled by this fake plateau into extracting a true ground state energy.It is required that the ground state energy is identified by much later time data t = O(10) fm in the temporal correlation functions, but this is very difficult due to the exponential growth of the noise over the signal.
We have introduced a consistency check (sanity check) to rule out obviously false results, combing the Lüscher's finite volume formula and the effective range expansion at low-energies.We check all the existing LQCD results obtained from the plateau method for the NN system, YKU2011, YIKU2012, YIKU2015, NPL2012, NPL2013, NPL2015 and CalLat2017.The results that none of them pass the test bring serious doubt on the existence of the bound states in the NN system claimed in those studies.In order to determine correct spectra of the NN system, at least, a much more sophisticated method than the simple plateau fitting, such as the variational method, must be employed.
We also have examined the reliability of the time-dependent HAL QCD method in which the contamination of the elastic excited states is well-tamed, unlike the plateau method.In time-dependent HAL QCD method, the source operator dependence can be handled in the velocity expansion, and the good convergence of the expansion is shown.We show the consistency between the Lüscher's finite volume formula and the HAL QCD method.Taking into account the finite volume energies obtained from the HAL QCD potential, it turns out that the above fake-plateau structure is indeed caused by the contamination of low-lying elastic excited states.
On the basis of the HAL QCD method, we have studied the πJ/ψ-ρη c -DD * coupled-channel interaction using (2+1)-flavor full QCD gauge configurations in order to reveal the nature of the tetraquark candidate Z c (3900).Thanks to the HAL QCD method, the full coupled-channel potential V αβ is obtained whose diagonal components are found to be small, so that Z c (3900) cannot be a simple πJ/ψ or DD * molecular state.Also, we find the strong off-diagonal couplings between πJ/ψ and DD * , which indicates that the Z c (3900) can be a threshold cusp.To confirm this, we calculate the invariant mass spectra and the pole position associated with the coupled-channel two-body S-matrix on the basis of V αβ .The results support that the peak in the πJ/ψ invariant mass spectrum is not associated with a conventional resonance state but is a threshold cusp induced by the strong πJ/ψ-DD * coupling.To further strengthen our conclusion, we perform a semi-phenomenological analysis of the three-body decay of the Y(4260), and find that the experimentally observed peak structures just above the DD * threshold are well reproduced in the Y(4260) → ππJ/ψ decay.

Figure 1 .
Figure 1.A schematic illustration of the finite volume spectrum in LQCD simulations for (a) single nucleon and (b) two nucleons.

Figure 2 .
Figure 2. The finite volume spectrum employed in the mock-up data.In the data, δE el = 50 MeV and δE inel = 500 MeV are taken, which are the typical lowest excited energies of the elastic and inelastic scattering states with L 4 fm and m π 500 MeV, respectively.
(c), we show the infinite volume extrapolations of the effective energies in the ΞΞ ( 3 S 1 ) channel using the data at L = 32, 40, 48, 64.The blue filled squares (red filled circles) are obtained from the smeared (wall) source, and the energy shifts ∆E ∞ at L → ∞ (open square and circle) are

Figure 4 .
Figure 4.The source operator dependence of the effective energy shift ∆E eff ΞΞ (t) in (a) ΞΞ ( 1 S 0 ) and (b) ΞΞ ( 3 S 1 ) channels at L = 48.The blue and red data denote ∆E eff (t) obtained from the smeared and wall sources, respectively.(c) The infinite volume extrapolation of the effective energies in the ΞΞ ( 3 S 1 ).The blue (red) data are obtained with the smeared (wall) source at L = 32, 40, 48, 64.(d) The sink operator dependence of the effective energy shift with the smeared source in the ΞΞ ( 1 S 0 ).The blue filled squares refer to the same data as in figure (a).The figures are taken from Ref. [30].

Figure 5 .
Figure 5.The schematic illustration for the relation between k cot δ(k)/m π and (k/m π ) 2 on finite volumes.The dashed lines are obtained from the Lüscher's formula for each volume.The discrete points are realized as satisfying both the Lüscher's formula and the ERE.The bound state pole condition is shown by filled circle.(a) The physical NN ( 3 S 1 ) system with empirical values of the ERE parameters.The ERE line is denoted by the red solid line.(b) An unphysical example of the relation between the Lüscher's formula and the ERE.The false results of ∆E lead to inconsistent ERE for positive (ERE k 2 >0, blue solid line) and for negative (ERE k 2 <0 , red solid line) energies.The effective range obtained from ERE k 2 <0 is also unnaturally large.These figures are provided by T. Iritani (See also Ref.[31]).
(a), we show the NN ( 3 S 1 ) ERE line at next-to-leading order (red solid line) and the bound state (pole) condition k cot δ(k) = − √ −k 2 (black solid line) with the empirical values of the scattering parameters, am π = −3.8 and r e m π = 1.3 for m π = 0.14 GeV.If the energy shift is reliably extracted on the finite box, k cot δ(k) obtained from Lüscher's formula (open triangles, squares and diamonds in figure 5 (a)) become consistent with the ERE line at low energies.Shown in figure 5 (b) is a typical example when one sees mirage of the temporal correlation function and thus misidentifies the plateau: the EREs at positive (ERE k 2 >0 , blue solid line) and negative (ERE k 2 <0 , red solid line) energies are inconsistent with each other, and the ERE k 2 <0

Figure 7 .
Figure 7. (a) The source operator dependence of the normalized correlation function given in Eq. (15) for time-slice t = 10-16.(b) The effective LO potential V eff calculated from R wall (r, t) in Eq. (19) (red) and the LO potential V LO obtained from both R smear (r, t) and R wall (r, t) in Eq. (20) (blue) at L = 64.

Figure 9 .
Figure 9. (a) The reconstructed effective energy shifts ∆E smear/wall eff (t) using three low-lying eigen-functions.(b) The contamination of the elastic excited states |b n /b 0 |.Filled (open) symbols represent the positive (negative) sign of b n /b 0 .The figures are taken from Ref. [33].
(b), in which the ratios |b n /b 0 | for both the smeared and wall sources are plotted.The values |b n /b 0 | describe the measure of the contamination of the n-th elastic excited states in the normalized correlation function R(r, t) defined in Eq. (
Figure 12.(a) The two-body invariant mass spectra in the πJ/ψ (red circles, scaled by 5), ρη c (green triangles) and DD * (blue squares) channels.The two-body πJ/ψ spectrum without the off-diagonal component of V αβ is also shown by Im f πJ/ψ,πJ/ψ 0 (black crosses, scaled by 25).(b) The pole of the S-matrix on the [bbb] sheet in the notation of [68] for πJ/ψ, ρη c and DD * channels (z = m α 1 + m α 2 + p 2 α /2µ α ).The connection to the physical [ttt] sheet is denoted by the red arrow (Re[z] ≥ m D + m D * ).Both figures correspond to case I in Table2.In figure12(a), the shorter error is statistical, while the longer one is statistical and systematic combined in quadrature.The figure is taken from Ref.[37].

Figure 14 .
Figure 14.The invariant mass spectrum of Y(4260) → ππJ/ψ calculated with V αβ for case I in table 2. The shaded area represents the statistical error.The vertical arrow shows the predicted peak position from the calculation.The gray dashed line shows the invariant mass spectrum of the Y(4260) decay without the off-diagonal components of V αβ .The figure is taken from Ref.[38].

Table 2 .
Meson masses calculated in the simulations in MeV unit.