# Charmoniumlike channels $1^{+}$ with isospin $1$ from lattice and effective
field theory

Mitja Sadl<EMAIL_ADDRESS>Faculty of Mathematics and Physics,
University of Ljubljana, Ljubljana, Slovenia  Sara Collins  Institut für
Theoretische Physik, Universität Regensburg, Regensburg, Germany  Zhi-Hui Guo
Department of Physics and Hebei Key Laboratory of Photophysics Research and
Application, Hebei Normal University, Shijiazhuang, China  M. Padmanath  The
Institute of Mathematical Sciences, CIT Campus, Chennai, 600113, India  Homi
Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai
400094, India  Sasa Prelovsek<EMAIL_ADDRESS>Faculty of Mathematics
and Physics, University of Ljubljana, Ljubljana, Slovenia  Jozef Stefan
Institute, Ljubljana, Slovenia  Lin-Wan Yan  Department of Physics and Hebei
Key Laboratory of Photophysics Research and Application, Hebei Normal
University, Shijiazhuang, China

###### Abstract

Many exotic charmoniumlike mesons have already been discovered experimentally,
of which the $Z_{c}$ mesons with isospin 1 are prominent examples. We
investigate $J^{PC}=1^{+\pm}$ states with flavor $\bar{c}c\bar{q}q$ ($q=u,d$)
in isospin 1 using lattice QCD. This is the first study of these mesons
employing more than one volume and involving frames with nonzero total
momentum. We utilize two $N_{f}=2+1$ CLS ensembles with
$m_{\pi}\simeq~{}$280\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. As the
simulations are performed with unphysical quark masses and at a single lattice
spacing of $a=$0.086\text{\,}\mathrm{f}\mathrm{m}$$, our results provide only
qualitative insights. Resulting eigenenergies are compatible or just slightly
shifted down with respect to noninteracting energies, where the most
significant shifts occur for certain $D\bar{D}^{*}$ states. Both channels
$1^{+\pm}$ have a virtual pole slightly below the threshold if $D\bar{D}^{*}$
is assumed to be decoupled from other channels. In addition, we perform a
coupled channel analysis of $J/\psi\pi$ and $D\bar{D}^{*}$ scattering with
$J^{PC}=1^{+-}$ within an effective field theory framework. The $J/\psi\pi$
and $D\bar{D}^{*}$ invariant-mass distributions from BESIII and finite-volume
energies from several lattice QCD simulations, including this work, are fitted
simultaneously. All fits yield two poles relatively close to the
$D\bar{D}^{*}$ threshold and reasonably reproduce the experimental $Z_{c}$
peaks. They also reproduce lattice energies up to slightly above the
$D\bar{D}^{*}$ threshold, while reproduction at even higher energies is better
for fits that put more weight on the lattice data. Our findings suggest that
the employed effective field theory can reasonably reconcile the peaks in the
experimental line shapes and the lattice energies, although those lie close to
noninteracting energies. We also study $J/\psi\pi$ scattering in s-wave and
place upper bounds on the phase shift.

##  I Introduction

The first glimpse into the realm of unconventional mesons happened with
Belle’s 2003 discovery of the $\chi_{c1}(3872)$ Choi  _et al._ (2003). While
its quantum numbers $I(J^{PC})=0(1^{++})$ initially align with a conventional
$\bar{c}c$ content, its mass and decay pattern suggest a more intricate
composition related to the nearby $D\bar{D}^{*}$ threshold. In 2013, the
BESIII and Belle collaborations unveiled the $Z_{c}^{+}(3900)$ with explicitly
exotic flavor composition $\bar{c}c\bar{d}u$ since it was observed in
$J/\psi\pi^{+}$ and $(D\bar{D}^{*})_{I=1}$ final states Ablikim  _et al._
(2013); Liu  _et al._ (2013a); Ablikim  _et al._ (2014a,
2015a)111$D\bar{D}^{*}$ generically refers to a linear combination of
$D\bar{D}^{*}$ and $\bar{D}D^{*}$ with appropriate charge parity $C=\pm$
throughout the text.. Recent determinations assert that the $Z_{c}(3900)$ is a
$1(1^{+-})$ state with a mass of $M=3887.1(2.6)\,\mathrm{MeV}$ and decay width
$\Gamma=28.4(2.6)\,\mathrm{MeV}$ Ablikim  _et al._ (2017a); Workman  _et al._
(2022). These charmoniumlike states reside near open-charm threshold
$D\bar{D}^{*}$, which is the energy region explored in this paper.

Various binding mechanisms have been proposed for the $Z_{c}(3900)$: it could
be a hadronic molecule Wang  _et al._ (2013); Dong  _et al._ (2013); Wilbring
_et al._ (2013); Aceti  _et al._ (2014); Gong  _et al._ (2016); Chen  _et al._
(2023); Liu  _et al._ (2024a, b), possess a compact tetraquark structure
Maiani  _et al._ (2013), or arise from a kinematic effect (either as threshold
cusp effect Liu and Li (2013); Swanson (2015, 2016) or triangle singularity
Liu  _et al._ (2016); Szczepaniak (2015); Chen  _et al._ (2023); Liu  _et al._
(2024a)). Although many articles propose it to be a molecular state, there are
no model-independent conclusions yet. Several binding mechanisms reasonably
reproduce $Z_{c}$ line shapes according to the analysis by JPAC Pilloni  _et
al._ (2017). Functional methods suggest it has significant $J/\psi\pi$ and
$D\bar{D}^{*}$ Fock components Hoffer  _et al._ (2024). A more extended
discussion of phenomenological studies is given in Sec. VIII.1.

Several lattice simulations of the $Z_{c}(3900)$ have been performed recently.
Two studies employing the HAL QCD approach Ikeda  _et al._ (2016); Ikeda and
for HAL QCD Collaboration (2018) highlighted the significance of cross-channel
interaction, aligning with the conclusions of Ref. Ortega  _et al._ (2019).
However, studies employing the Lüscher formalism could not confirm a (narrow)
resonance-like peak near the threshold. This includes Prelovsek and Leskovec
(2013a); Prelovsek  _et al._ (2015); Chen  _et al._ (2014); Lee  _et al._
(2014); Cheung  _et al._ (2017), as well as the coupled channel analysis of
Chen  _et al._ (2019). Notably, no additional finite-volume eigenstates were
found and the energy shifts relative to the noninteracting levels were
insignificant. Direct comparisons between methods are challenging as the HAL
QCD approach does not provide information on the energy shifts.

While experimental observations have revealed charmoniumlike states with
$I(J^{PC})=1(1^{+-})$, no $\bar{c}c\bar{q}q$ states with $1(1^{++})$ have been
detected.222In the following,
$\bar{q}q\in\\{\bar{u}d,\bar{u}u-\bar{d}d,\bar{d}u\\}$ since $I=1$ is
considered. Such a state would be an isospin partner of $\chi_{c1}(3872)$. An
isospin-1 partner was predicted by diquark-antidiquark models Maiani  _et al._
(2005); Anwar  _et al._ (2018), and there are arguments for its existence also
in the molecular scenario Ji  _et al._ (2022); Zhang  _et al._ (2024), where
reasons for non-observation in experiment are provided Zhang  _et al._ (2024).
Two lattice QCD studies Prelovsek and Leskovec (2013b); Padmanath  _et al._
(2015), which found the state $\chi_{c1}(3872)$ slightly below the
$D\bar{D}^{*}$ threshold, also did not find significant interaction in the
isospin-1 system. However, the scattering amplitude for the $1(1^{++})$
channel has not yet been extracted from the lattice, and this is one of the
aims of the present paper.

This article presents the results of a lattice study of channels with the
flavor composition $\bar{c}c\bar{q}q$ and the quantum numbers
$I(J^{PC})=1(1^{+\pm})$.333In practice, we implement the charge-neutral
($I_{z}=0$) combination $\bar{u}u-\bar{d}d$. Our central result is the
extraction of the discrete spectrum of QCD eigenstates within a finite volume.
To realize this, we compute matrices of two-point correlation functions using
a large basis of interpolators. We construct the interpolators from pairs of
meson operators (each meson operator with definite momentum) that are
projected onto two different total momenta,
$|{\hbox{\boldmath{$P$}}}|L/(2\pi)=0,1$. The matrices of correlation functions
are evaluated on two ensembles with different spatial extents and numerous
excited finite-volume energies are extracted. In the $1^{++}$ channel, the
$J/\psi\rho$ threshold is positioned below the $D\bar{D}^{*}$ threshold. The
$1^{+-}$ channel features two thresholds, $J/\psi\pi$ and $\eta_{c}\rho$,
below $D\bar{D}^{*}$. These lower-lying thresholds result in inelastic
$D\bar{D}^{*}$ scattering and introduce significant complexity to extracting
the pertinent $D\bar{D}^{*}$ eigenenergies and scattering matrices.

Resonances and near-threshold (virtual) bound states manifest as peaks in
hadron scattering amplitudes. These amplitudes are computed on the lattice
(and also extracted from experimental measurements) at real values of the
scattering energy. A comprehensive understanding of these amplitudes comes
from considering their singularity structure in the complex energy plane.
Scattering amplitudes can exhibit pole singularities that may be identified
with bound states and virtual states (situated on the real axis below the
lowest threshold444Scattering particles have positive (negative) imaginary
momenta at the energy of the bound (virtual) state.) or resonances (positioned
off the real axis).

After obtaining finite-volume energies from two distinct volumes and two
frames, we attempt to determine scattering amplitudes. Extracting the energy
dependence of these amplitudes is intricate. In this work, we perform the
scattering analysis following two different approaches, as detailed in the
following two paragraphs.

In the first approach, we try to address the question of whether a
charmoniumlike, isospin 1 four-quark state with $1^{++}$ or $1^{+-}$ is
present near the $D\bar{D}^{*}$ threshold with our lattice setup. We aim to
discern whether this state exists, particularly in the scenario where the
interaction between $D\bar{D}^{*}$ and other channels has a negligible impact.
Therefore, scattering amplitudes are fitted from our eigenenergies using the
Lüscher method assuming one-channel $D\bar{D}^{*}$ scattering. The energy
dependence of scattering amplitude is parametrized using an effective range
expansion.

The second approach is applied only to $1^{+-}$: the
$J/\psi\pi$-$D\bar{D}^{*}$ coupled-channel system within a covariant framework
of a contact effective field theory (EFT) is studied. More precisely, contact
interactions $D\bar{D}^{*}$-$D\bar{D}^{*}$ and $D\bar{D}^{*}$-$J/\psi\pi$ are
considered. We perform a joint fit to the experimental $J/\psi\pi$ and
$D\bar{D}^{*}$ invariant-mass distributions from BESIII Ablikim  _et al._
(2017a, 2015b), our and other Cheung  _et al._ (2017); Chen  _et al._ (2019)
lattice finite-volume energy levels. The resulting coupled-channel scattering
amplitude features resonance poles near $D\bar{D}^{*}$ threshold. A similar
investigation was carried out in Ref. Yan  _et al._ (2024), where the lattice
results of this work were not included.

This paper is structured as follows. The next three sections discuss the
lattice setup, the creation/annihilation operators and the extraction of
finite-volume energy levels. These levels are presented in Sec. V. Scattering
amplitudes in the one-channel approximation are extracted in Sec. VI. Finite-
volume energies are fitted together with the experimental data within a
covariant framework of a coupled-channel system in Sec. VII. Other theoretical
studies are discussed in Sec. VIII. The practically vanishing $J/\psi\pi$
energy shifts indicate a negligible interaction between $J/\psi$ and $\pi$,
and we put the bound on the $J/\psi\pi$ scattering length in Sec. IX. This is
followed by the conclusions.

##  II Lattice details

We employ two ensembles of gauge field configurations generated with
$N_{\textrm{f}}=2+1$ non-perturbatively $\mathcal{O}(a)$ improved Wilson
dynamical fermions, provided by the Coordinated Lattice Simulations (CLS)
consortium Bruno  _et al._ (2015); Bali  _et al._ (2016). The quark masses
$m_{u/d,s}$ are chosen to lie on a trajectory that approaches the physical
point holding the flavor average quark mass, $2m_{u/d}+m_{s}$, constant. The
ensembles have volumes $N_{T}\times N_{L}^{3}=128\times 24^{3}$ (denoted U101)
and $96\times 32^{3}$ (labeled H105). We utilize $255$ and $492$
configurations on the smaller and larger volumes, respectively Bruno  _et al._
(2017). Both ensembles have an inverse gauge coupling
$\beta=6/g^{2}=$3.4\text{\,}$$, lattice spacing $a=0.08636(98)(40)\,$fm and
pion mass $m_{\pi}=$280(3)\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. Open
boundary conditions in time are imposed Luscher and Schaefer (2013) and the
sources and sink timeslices of the correlation functions are located in the
bulk away from the boundary.

The study is performed with a slightly larger than physical charm quark mass
Piemonte  _et al._ (2019). The masses of the pion, $D$, $D^{*}$ mesons and the
spin-averaged 1S-charmonium mass determined on the larger ensemble are shown
in Table 1. The hadron masses on both ensembles and the corresponding
experimental values are given in Appendix A. These masses indicate that the
$D^{*}$ meson is stable on our lattices.

Table 1: Masses of the pion, $D$, $D^{*}$ mesons and $M_{\textrm{av}}=\left(m_{\eta_{c}}+3m_{J/\psi}\right)/4$ measured on the larger lattice with $N_{L}=32$. $m_{\pi}$ [MeV] |  $m_{D}$ [MeV] |  $m_{D^{*}}$ [MeV] |  $M_{\textrm{av}}$ [MeV]
---|---|---|---
$280(3)$ |  $1927(1)$ |  $2049(2)$ |  $3104.3(4)$
Table 2: Our setup for the $\bar{c}c\bar{q}q$ system with total momenta ($P$), spatial lattice symmetry groups (LG), lattice irreps ($\Lambda^{(P)}$), total angular momentum–parities ($J^{P}$) and partial waves ($\ell$) that contribute to each $\Lambda^{(P)}$. For every combination of $\Lambda^{(P)}$, charge parity ($C$) and the number of lattice sites in spatial direction ($N_{L}$), the total number of interpolators used is presented. We list only $J\leq 3$ and $\ell\leq 2$. We do not calculate the correlation function for the case $\Lambda=A_{2}$ with $C=-$ and larger $N_{L}=32$ since it carries the most populated spectrum, and extraction of its energy levels would be too noisy. $P$ | LG | $\Lambda^{(P)}$ | J^P |  $\ell$ | $C$ | (less) important interpolator types | $N_{L}$ | total # of interpolators
---|---|---|---|---|---|---|---|---
$(0,0,0)$ | O${}_{\textrm{h}}$ | $T_{1}^{+}$ | 1^+,3^+ |  $0,2$ | $-$ | $J/\psi\pi,\eta_{c}\rho,D\bar{D}^{*}(D^{*}\bar{D}^{*},h_{c}\pi)$ | 24 | 15
|  |  |  |  |  |  | 32 | 21
|  |  |  |  | $+$ | $J/\psi\rho,D\bar{D}^{*}(\chi_{c0}\pi,\chi_{c1}\pi)$ | 24 | 5
|  |  |  |  |  |  | 32 | 10
$(0,0,1)\frac{2\pi}{L}$ | Dic4 | $A_{2}$ |  |  $0,1,2$ | $-$ | $J/\psi\pi,\eta_{c}\rho,D\bar{D}^{*}$ | 24 | 21
|  |  | 0^-,1^+, |  |  |  | 32 | /
|  |  | 2^-,3^+ |  | $+$ | $J/\psi\rho,D\bar{D}^{*}(\chi_{c0}\pi,\chi_{c1}\pi,\eta_{c}a_{0})$ | 24 | 13
|  |  |  |  |  |  | 32 | 17

The Wick contractions are evaluated using the full Peardon  _et al._ (2009)
and stochastic Morningstar  _et al._ (2011) distillation method with 90 and
100 Laplacian eigenvectors for $N_{L}=24$ and $N_{L}=32$, respectively. The
setup for calculating the quark propagators is the same as that employed in
Refs. Piemonte  _et al._ (2019); Prelovsek  _et al._ (2021), except for the
number of eigenvectors utilized. The correlation matrices are averaged over
all three polarizations (momentum directions) and up to nine source-time
slices to increase the statistical precision. All uncertainties quoted are
statistical only and are determined using the bootstrap method. We quote
symmetric and asymmetric uncertainties, as defined in appendix A of Ref.
Prelovsek  _et al._ (2021), using $N_{\textrm{b}}=1000$ bootstrap samples.
Since this investigation utilizes lattice gauge ensembles featuring two
distinct physical volumes, a single lattice spacing, and nonphysical quark
masses, it is essential to note that our findings allow only a qualitative
comparison with the experimental data.

##  III Two-hadron interpolating operators with $I=1$

The standard approach for determining finite-volume energies is employed,
wherein we calculate two-point correlation functions of the form

$C_{ij}(t)=\langle 0|O_{i}(t)O_{j}^{\dagger}(0)|0\rangle\>.$ (1)

Here, $O_{i}$ ($O_{j}^{\dagger}$) is an operator that annihilates (creates) a
state with specific quantum numbers. These interpolating operators555We
interchangeably use notions operator and interpolator for simplicity. are
constructed as a product of one or more quark and antiquark fields. We
implement only meson-meson type interpolators

$\bar{c}\Gamma_{i}c(|{\hbox{\boldmath{$p$}}}_{i}|^{2})\
\bar{q}\Gamma_{j}q(|{\hbox{\boldmath{$p$}}}_{j}|^{2})\quad\textrm{and}\quad\bar{c}\Gamma_{i}q(|{\hbox{\boldmath{$p$}}}_{i}|^{2})\
\bar{q}\Gamma_{j}c(|{\hbox{\boldmath{$p$}}}_{j}|^{2})\>.$ (2)

One-meson interpolators
$\sum_{{\hbox{\boldmath{$x$}}}}e^{i{\hbox{\boldmath{$p$}}}\cdot{\hbox{\boldmath{$x$}}}}\bar{\psi}({\hbox{\boldmath{$x$}}})\Gamma\psi({\hbox{\boldmath{$x$}}})$
have an appropriate Dirac structure $\Gamma$ and are separately projected to a
definite momenta, where the total momentum is
${\hbox{\boldmath{$P$}}}={\hbox{\boldmath{$p$}}}_{i}+{\hbox{\boldmath{$p$}}}_{j}$.666Throughout
the article, ${\hbox{\boldmath{$p$}}}_{i,j}$ and
$\left|{\hbox{\boldmath{$p$}}}_{i,j}\right|^{2}$, which represent the momentum
of a single meson, will be given in units of $2\pi/L$ and $(2\pi/L)^{2}$,
respectively, since
${\hbox{\boldmath{$p$}}}_{i,j}={\hbox{\boldmath{$n$}}}_{i,j}\frac{2\pi}{L}$,
with ${\hbox{\boldmath{$n$}}}_{i,j}\in N_{L}^{3}$. One-hadron interpolators of
type $\bar{c}\Gamma c$ are absent in the $I=1$ sector. We omit local diquark-
antidiquark interpolators since previous lattice studies Cheung  _et al._
(2017); Padmanath  _et al._ (2015); Prelovsek  _et al._ (2015) suggest they
have little influence on the rest of the spectrum. Some meson-meson channels
represented by the left interpolator in (2) include effects from light
resonances, specifically the $\rho$ and $a_{0}$, which are unstable at our
pion mass. Despite this, we only implement two-hadron interpolators; no three-
hadron interpolators are included. The latter would be more suitable for
capturing the physics of decays involving three hadrons in the final state,
e.g., $Z_{c}(3900)\to\eta_{c}\rho\to\eta_{c}\pi\pi$.

For $I=1$, the Wick contractions for the matrix of correlation functions (Eq.
(1)) involve diagrams where the light quarks exclusively propagate from the
source to the sink. Two types of quark line contractions arise for the charm
quarks: one where the charm quarks propagate from the source to the sink
timeslice and another where they propagate from and to the same source/sink
timeslice. The latter annihilation diagrams can mix with channels that do not
contain charm quarks. Similar to most previous works, we omit these diagrams
as it is very challenging to include this mixing in the lattice calculation.
We remark that experiments do not observe these decay channels in the region
of interest.

The setup for all the interpolators implemented in this study is provided in
Table 2. Their free energies, in the center-of-momentum (cm) frame, mainly
span the region below $4.2\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. For more
comprehensive lists, we refer to Tables 7 and 8 in Appendix B. To study
$s$-wave $D\bar{D}^{*}$ scattering with spin-parity $1^{+}$, we focus on two
finite-volume irreducible representations (irreps) of the spatial lattice
symmetry groups $O_{h}$ ($|{\hbox{\boldmath{$P$}}}|=0$) and $\textrm{Dic}_{4}$
($|{\hbox{\boldmath{$P$}}}|=1\cdot 2\pi/L$), namely $\Lambda^{P}=T_{1}^{+}$
and $\Lambda=A_{2}$, respectively.

Figure 1: Eigenenergies in the rest frame ($|{\hbox{\boldmath{$P$}}}|=0$)
finite-volume irrep $T_{1}^{+-}$ relevant for the quantum numbers
$I(J^{PC})=1(1^{+-})$. The left and right panes represent the energy levels
(points, $E_{n,\textrm{cm}}$) and noninteracting energies (lines,
$E_{\mathrm{ni,cm}}^{\mathrm{con}}$) for the case where all interpolators are
included and the $\eta_{c}\rho$ interpolators are excluded, respectively.
States that dominantly couple to $D\bar{D}^{*}$ ($J/\psi\pi$) interpolators
are presented with relatively warm (cool) colors. Empty gray symbols and
dashed gray lines represent states with dominant overlap to $\eta_{c}\rho$
interpolators containing the $\rho$ resonance, which in our case is treated as
stable under the strong interaction. Numbers within the square brackets refer
to the multiplicity of certain noninteracting levels. If there is a
possibility to relate an extracted energy level to a certain partial wave,
then its symbol is a square, diamond or triangle for s-wave ($J=1$), d-wave
($J=1$) or d-wave ($J=3$), respectively. Finite-volume eigenenergies are shown
with $1\sigma$ statistical uncertainty, and some are slightly horizontally
shifted for clarity. The main text explains that the noninteracting lines are
drawn with each meson mass being a linear function of the lattice size,
$m(N_{L})$.

Figure 2: Eigenenergies in the moving frame
($|{\hbox{\boldmath{$P$}}}|L/(2\pi)=1$) finite-volume irrep $A_{2}^{-}$
relevant for the quantum numbers $I(J^{PC})=1(1^{+-})$. The left and right
panes represent the energy levels (points, $E_{n,\textrm{cm}}$) and
noninteracting energies (lines, $E_{\textrm{ni,cm}}^{\textrm{con}}$) for the
case where all interpolators are included and the $\eta_{c}\rho$ interpolators
are excluded, respectively. Our calculations in this irrep are performed only
on the smaller lattice (with $N_{L}=24$). For more information see the caption
of Fig. 1

Figure 3: Eigenenergies in the rest frame ($|{\hbox{\boldmath{$P$}}}|=0$)
finite-volume irrep $T_{1}^{++}$ relevant for the quantum numbers
$I(J^{PC})=1(1^{++})$. The left and right panes represent the energy levels
(points, $E_{n,\textrm{cm}}$) and noninteracting energies (lines,
$E_{\textrm{ni,cm}}^{\textrm{con}}$) for the case where all interpolators are
included and the $J/\psi\rho$ and $\eta_{c}a_{0}$ interpolators are excluded,
respectively. States that dominantly couple to $D\bar{D}^{*}$ interpolators
(here with $C=+$ and not $C=-$) are presented with warm colors. Empty gray
symbols and dashed gray lines represent states coupled to $J/\psi\rho$
interpolators containing the $\rho$ resonance, which is treated as stable in
our case. For more information see caption of Fig. 1.

Figure 4: Eigenenergies in the moving frame
($|{\hbox{\boldmath{$P$}}}|L/(2\pi)=1$) finite-volume irrep $A_{2}^{+}$
relevant for the quantum numbers $I(J^{PC})=1(1^{++})$. The left and right
panes represent the energy levels (points, $E_{n,\textrm{cm}}$) and
noninteracting energies (lines, $E_{\textrm{ni,cm}}^{\textrm{con}}$) for the
case where all interpolators are included and the $J/\psi\rho$ and
$\eta_{c}a_{0}$ interpolators are excluded, respectively. Empty symbols and
dashed lines represent states coupled to $J/\psi\rho$ and $\eta_{c}a_{0}$
interpolators containing the $\rho$ and $a_{0}$ resonances, which in our case
are treated as stable. As commented in the main text, we present the
$\eta_{c}a_{0}$ noninteracting energies discontinuously since the “masses” of
$a_{0}$ measured on both ensembles differ significantly. For more information
see the caption of Fig. 1. Figure 5: A close-up of the densest region of the
spectrum for irrep $A_{2}^{+}$ on the $N_{L}=32$ ensemble shown in the left-
hand side plot in Fig. 4. Points are shifted horizontally for clarity.

The partial-wave method and the projection method (see Ref. Prelovsek  _et
al._ (2017)) are employed for the construction of interpolators in the case of
$|{\hbox{\boldmath{$P$}}}|=0$ and $|{\hbox{\boldmath{$P$}}}|=1\cdot 2\pi/L$,
respectively. In the absence of interaction, two-particle states can exhibit
energy degeneracy if at least one particle has a non-zero spin. This
degeneracy depends on the spin and momentum of both particles. As an
illustrative example, consider the case of $J/\psi\pi$ (vector-pseudoscalar)
interpolators that transform according to a certain “row” of the 3D irrep
$T_{1}^{+-}$. The total spin of a pseudoscalar-vector two-meson system is
$S=1$. Considering $J^{P}=1^{+},\ 3^{+},\ \ldots\ $ that contribute to
$T_{1}^{+-}$, one can have the partial waves $\ell=0,\ 2,\ \ldots\ $. When
both particles are at rest, the orbital angular momentum $\ell=0$ (s-wave)

$\displaystyle O^{J/\psi(0)\pi(0)}_{S=1,\ell=0,J=1}$
$\displaystyle=\bar{c}\gamma_{z}c({\hbox{\boldmath{$0$}}})\
\bar{q}\gamma_{5}q({\hbox{\boldmath{$0$}}})\>.$ (3)

Considering $J/\psi(1)\pi(1)$, the $\ell=0$ as well as $\ell=2$ (d-wave)
appear and the partial-wave projected interpolators are Prelovsek  _et al._
(2017)

$\displaystyle O^{J/\psi(1)\pi(1)}_{S=1,\ell=0,J=1}$
$\displaystyle=\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{x})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{x})\,+\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{x})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{x})$
$\displaystyle+\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{y})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{y})\,+\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{y})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{y})$
$\displaystyle+\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{z})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{z})\,+\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{z})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{z})$
$\displaystyle O^{J/\psi(1)\pi(1)}_{S=1,\ell=2,J=1}$
$\displaystyle=\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{x})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{x})\,+\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{x})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{x})$
$\displaystyle+\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{y})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{y})\,+\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{y})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{y})$
$\displaystyle-2\,\bar{c}\gamma_{z}c(\hat{{\hbox{\boldmath{$e$}}}}_{z})\,\bar{q}\gamma_{5}q(-\hat{{\hbox{\boldmath{$e$}}}}_{z})$
$\displaystyle-2\,\bar{c}\gamma_{z}c(-\hat{{\hbox{\boldmath{$e$}}}}_{z})\,\bar{q}\gamma_{5}q(\hat{{\hbox{\boldmath{$e$}}}}_{z})\>.$
(4)

Most of our interpolators consist of pseudoscalar (P) and vector (V) mesons.
As shown in Eq. (4), the case $P(1)V(1)$ exhibits a multiplicity of 2 linearly
independent interpolators, $O^{P(1)V(1)}_{S=1,\ell=0,J=1}$ and
$O^{P(1)V(1)}_{S=1,\ell=2,J=1}$. A similar degeneracy appears when considering
$P(3)V(3)$ interpolators, $O^{P(3)V(3)}_{S=1,\ell=0,J=1}$ and
$O^{P(3)V(3)}_{S=1,\ell=2,J=1}$. Furthermore, there are 3 linearly independent
$P(2)V(2)$ interpolators: $O^{P(2)V(2)}_{S=1,\ell=0,J=1}$,
$O^{P(2)V(2)}_{S=1,\ell=2,J=1}$ and $O^{P(2)V(2)}_{S=1,\ell=2,J=3}$. The
$P(1)V(2)$, $P(2)V(1)$ and $V(1)V(1)$ energy levels have degeneracies of 2, 2
and 3, respectively.

Our objective is to extract all these eigenstates, including those that are
degenerate in the noninteracting limit. Therefore, for every given
interpolator type, we implement all its linearly independent combinations.

##  IV Extraction of finite-volume energy levels

On a periodic lattice with spatial size $L=N_{L}a$, a two-hadron system
exhibits discrete energies in noninteracting (ni) limit

$E_{\textrm{ni}}^{\textrm{lat}}=E_{H_{i}({\hbox{\boldmath{$p$}}}_{i})}^{\textrm{lat}}+E_{H_{j}({\hbox{\boldmath{$p$}}}_{j})}^{\textrm{lat}}\>.$
(5)

In the continuum limit, the noninteracting energies follow the relativistic
dispersion relation, and

$E_{\textrm{ni}}^{\textrm{con}}=E_{H_{i}({\hbox{\boldmath{$p$}}}_{i})}^{\textrm{con}}+E_{H_{j}({\hbox{\boldmath{$p$}}}_{j})}^{\textrm{con}}\>,\
E_{H_{i}({\hbox{\boldmath{$p$}}}_{i})}^{\textrm{con}}=\sqrt{m^{2}_{H_{i}}+\left|{\hbox{\boldmath{$p$}}}_{i}\right|^{2}}\>.$
(6)

Throughout this article, these energies are presented in the cm-frame,

$E_{\textrm{ni,cm}}^{\textrm{con}}=\sqrt{\left(E_{\textrm{ni}}^{\textrm{con}}\right)^{2}-\left|{\hbox{\boldmath{$P$}}}\right|^{2}}\>.$
(7)

In the interacting theory, finite-volume eigenenergies $E_{n}^{\textrm{lat}}$
are extracted from ab-initio lattice simulations. Specifically, these energies
are obtained through single-exponential fits to the eigenvalues
$\lambda^{(n)}(t,t_{0})\propto e^{-E_{n}^{\textrm{lat}}(t-t_{0})}$ of the
generalized eigenvalue problem (GEVP)

$C(t){\hbox{\boldmath{$u$}}}^{(n)}(t)=\lambda^{(n)}(t,t_{0})C(t_{0}){\hbox{\boldmath{$u$}}}^{(n)}(t)\>,$
(8)

where $t_{0}/a\in\\{2,3\\}$ have been used Michael (1985). Additional details
about our specific implementation can be found in Appendices C and D.

We emphasize that the finite-volume spectra are determined at a single lattice
spacing. Consequently, we cannot quantify the uncertainty arising from lattice
discretization effects, which may lead to the energy levels differing from
continuum expectations. These lattice artifacts are more significant for charm
quarks than light quarks since $am_{c}$ is not small for our lattices. The
measured single-hadron energies
$E_{H_{i}({\hbox{\boldmath{$p$}}}_{i})}^{\textrm{lat}}$ deviate slightly from
their continuum counterparts
$E_{H_{i}({\hbox{\boldmath{$p$}}}_{i})}^{\textrm{con}}$, especially for the
smaller volume where the unit of lattice momentum is larger. To mitigate these
effects, we first evaluate the energy shift of each interacting eigenstate
with respect to the noninteracting state

$\Delta E_{n}=E_{n}^{\textrm{lat}}-E_{\textrm{ni}}^{\textrm{lat}}\>.$ (9)

Subsequently, we determine corrected finite-volume eigenenergies using Eq. (6)

$E_{n}=\Delta E_{n}+E_{\textrm{ni}}^{\textrm{con}}$ (10)

and present them in the cm-frame

$E_{n,\textrm{cm}}=\sqrt{E_{n}^{2}-\left|{\hbox{\boldmath{$P$}}}\right|^{2}}\>.$
(11)

Note that the energies $E_{n}$ are, by construction, equal to
$E_{n}^{\textrm{lat}}$ in the continuum limit ($a\rightarrow 0$).

This procedure for extracting $E_{n}$ in Eq. (10) requires a clear
determination of the relevant noninteracting level for the eigenstate $n$. A
possible incorrect identification could modify $E_{n}$; however, the resulting
differences typically remain within the presented uncertainties, as
demonstrated in an example in Appendix E.

##  V Main results: finite-volume eigenenergies

The interacting energy levels evaluated using Eq. (11) are depicted in Figs.
1, 2 for $J^{PC}=1^{+-}$ and Figs. 3, 4 for $1^{++}$. These plots also include
the noninteracting energies determined using Eq. (7) presented as lines.
Numerous eigenstates feature in our spectra, which become increasingly dense
for non-zero total momentum and for the larger spatial volume. In particular,
many energy levels reside in the region above the
$D(1)\bar{D}^{*}(0)$/$D(0)\bar{D}^{*}(1)$ energies in the irrep $A_{2}^{+}$ on
the $N_{L}=32$ lattice and, for clarity, the spectrum for this reduced energy
region is presented in Fig. 5.

We use numbers inside square brackets in the figures to indicate the
degeneracy of energy levels in the noninteracting limit. The origin of this
degeneracy is discussed in Sec. III. Each linearly independent interpolator
yields a corresponding finite-volume interacting eigenenergy. Furthermore, the
states we extract typically couple well to interpolators constructed using the
partial-wave method.

Although the light mesons $\rho$ and $a_{0}$ are resonances, in our analysis,
we treat them as stable under the strong interaction. Note that a reliable
extraction of their eigenenergies would necessitate the inclusion of
3-particle operators, which is beyond the scope of the present work. As a
result, our analysis does not provide complete spectra, but our findings offer
insights into whether the energies of these “light mesons” are modified in the
presence of charmonium. In Figs. 1–4 (as indicated by the empty symbols), we
observe significant uncertainties in the eigenenergies associated with
interpolators containing $\rho$ and $a_{0}$.777We refrain from representing
the noninteracting energies related to the $a_{0}$ meson as continuous lines.
A significant difference exists between $m_{a_{0}}$ extracted on the
$N_{L}=24$ and $N_{L}=32$ ensembles. The unexpectedly low value for
$m_{a_{0}}$ on the smaller volume likely represents the sum of the masses of
the $\pi$ and $\eta$ mesons, which are allowed strong decay products of the
$a_{0}$.

The lines representing the noninteracting energies where both mesons are at
rest would appear horizontal if the meson masses were independent of $L$. This
is not valid exactly in our simulation; therefore, the lines for
noninteracting energies are interpolated linearly between the $N_{L}=24$ and
$32$ values in the figures.

In the $J^{PC}=1^{+-}$ sector, states coupled to $J/\psi\pi$ (apart from the
highest levels) align with the noninteracting expectations, suggesting that
interaction between $J/\psi$ and $\pi$ is very small. This topic is covered in
Sec. IX, where an upper bound on the scattering length of elastic s-wave
$J/\psi\pi$ scattering is estimated.

Our main objective is to find out whether there is any attraction888An
attraction manifests as an energy level lying significantly below the free
energy. between $D$ and $\bar{D}^{*}$ for both the $1^{+-}$ and $1^{++}$
channels, specifically when employing the full correlation matrix (as shown in
the left-hand side plots of Figures 1, 2, 3, and 4). The states that couple
predominantly to $D(0)\bar{D}^{*}(0)$, s-wave $D(1)\bar{D}^{*}(1)$,
$D(0)\bar{D}^{*}(1)$ and $D(1)\bar{D}^{*}(0)$ exhibit slightly negative energy
shifts. Exceptions are the states in the $1^{+-}$ channel on the $N_{L}=24$
ensemble in the $T_{1}^{+-}$ irrep that couple strongly to
$D(0)\bar{D}^{*}(0)$ and s-wave $D(1)\bar{D}^{*}(1)$, as well as $A_{2}^{-}$
levels that couple to $D(1)\bar{D}^{*}(0)$.

To qualitatively check how strongly the $J/\psi\pi$ and $D\bar{D}^{*}$
channels couple to $\eta_{c}\rho$ for $1^{+-}$, we omit the latter set of
interpolators and investigate the differences in the resulting spectrum.
Energy levels are extracted from correlation matrices that are submatrices of
the original ones and do not contain $\eta_{c}\rho$ interpolators. These
levels are presented on the right-hand side in Figs. 1, 2 and the
corresponding $\eta_{c}\rho$ energy levels do not feature in these spectra.
Apart from that, the spectra from the left- and right-hand side plots agree
within uncertainties, except for the s-wave $D(1)\bar{D}^{*}(1)$ state (for
the $N_{L}=32$ ensemble), which exhibits a greater degree of attraction when
the $\eta_{c}\rho$ interpolators are included.

A similar procedure is carried out for $1^{++}$, where eigenenergies are
extracted from submatrices that do not contain $J/\psi\rho$ and
$\eta_{c}a_{0}$ interpolators, whereas the $D\bar{D}^{*}$ interpolators are
still present. To ascertain the impact of the aforementioned interpolators,
compare the left- and right-hand sides in Figs. 3, 4. The conclusion is
similar to that for the $1^{+-}$ sector; the spectra agree within
uncertainties.

##  VI Approximation of one-channel s-wave $D\bar{D}^{*}$ scattering

Figure 6: $D\bar{D}^{*}$ scattering in partial wave $\ell=0$ in the one-
channel approximation with $J^{PC}=1^{+-}$. All interpolators were included in
the variational analysis on the left, while the $\eta_{c}\rho$ interpolators
were excluded on the right. Finite-volume energy levels incorporated into the
left and right fit (see the legend) can be found in Figs. 1, 2 in the left and
right plots, respectively. The black solid line displays $ka\cot{\delta_{0}}$
(see Eq. (14)) evaluated using the best fit parameters. The magenta octagon at
the intersection of the black and dashed cyan line (where $ik=+|k|$)
represents the presence of a virtual state pole. The gray vertical dashed
lines indicate the position of the branch point of the left-hand cut.

Figure 7: $D\bar{D}^{*}$ scattering in partial wave $\ell=0$ in the one-
channel approximation with $J^{PC}=1^{++}$. All interpolators were included in
the variational analysis on the left, while the $J/\psi\rho$ and
$\eta_{c}a_{0}$ were excluded on the right. Finite-volume energy levels
incorporated into the left and right fit (see the legend) can be found in
Figs. 3, 4 in the left and right plots, respectively. For more information see
the caption of Fig. 6.

We aim to calculate hadron scattering amplitudes using lattice QCD and examine
their singularity content. The former can be done using the Lüscher finite-
volume formalism and its extensions Lüscher (1986); Luscher (1991); Lüscher
(1991); Briceño (2014) (see Ref. Briceño  _et al._ (2018) for a review), which
relate finite-volume energies with the infinite-volume scattering amplitudes.
In this section, we assume that the $D\bar{D}^{*}$ channel is decoupled and
concentrate on single-channel $D\bar{D}^{*}$ scattering in partial wave
$\ell=0$ near the threshold for both $J^{PC}=1^{+-}$ and $1^{++}$. Our goal is
to check whether a charmoniumlike isospin 1 four-quark state with
$J^{P}=1^{+}$ close to $D\bar{D}^{*}$ threshold exists if the coupling of
$D\bar{D}^{*}$ to other channels is negligible. Note that this strategy cannot
determine whether the cross-channel interactions substantially affect an
exotic four-quark state. Possible effects from the left-hand cut (lhc) arising
from one-pion exchange are omitted in our analysis, as detailed at the end of
this section.

Figure 8: The quantity $k\left|t_{0}\right|^{2}$ proportional to the rate for
$D\bar{D}^{*}$ scattering in the one-channel approximation for $J^{PC}=1^{+-}$
(left) and $1^{++}$ (right). Here, all interpolators were included in the
variational analysis. The black line represents the central value
corresponding to the effective range parameters in Table 3. The gray band
shows the significant 1$\sigma$ uncertainty arising from the uncertainty of
the virtual pole position presented in violet.

For partial wave $\ell$, the scattering amplitude $t_{\ell}$ is defined
through the $S$-matrix,
$S_{\ell}=e^{2i\delta_{\ell}}=1+i\frac{4k}{E_{\textrm{cm}}}t_{\ell}$, where

$t_{\ell}=\frac{E_{\textrm{cm}}}{2}\frac{1}{k\cot{\delta_{\ell}}-ik},$ (12)

$\delta_{l}$ is the phase shift and $k$ is the magnitude of the spatial
momentum in the cm-frame, derived from
$E_{\textrm{cm}}=\sqrt{m_{D}^{2}+k^{2}}+\sqrt{m_{D^{*}}^{2}+k^{2}}$. The total
angular momentum is $J=1$, and the partial wave is $\ell=0$.999All energy
levels significantly overlapping with d-wave $D\bar{D}^{*}$ operators align
consistently with their corresponding noninteracting energies. Contributions
from partial waves with $\ell>0$ near the $D\bar{D}^{*}$ threshold—our region
of interest—are effectively suppressed by the phase space factor $k^{2\ell}$.
Therefore, we assume $t_{\ell\geq 2}=0$ and consider the mixing between
$\ell=2$ and $\ell=0$ in $J=1$ negligible. Utilizing Lüscher relation Briceño
(2014), we evaluate

$k\cot{\delta_{0}(k)}=\frac{2\mathcal{Z}^{{\hbox{\boldmath{$d$}}}}_{00}\left(1;(k\frac{L}{2\pi})^{2}\right)}{\pi^{1/2}\gamma
L}\>,$ (13)

where $\mathcal{Z}^{{\hbox{\boldmath{$d$}}}}_{00}$ is the Lüscher zeta
function Luscher (1991), and
${\hbox{\boldmath{$d$}}}={\hbox{\boldmath{$P$}}}\frac{L}{2\pi}$ is the
normalized total momentum. Finally, the following parametrization is chosen
for energies near the threshold:

$k\cot{\delta_{0}}=\frac{1}{a_{0}}+\frac{r_{0}}{2}k^{2}\>,$ (14)

which is an effective range expansion in $k^{2}$ up to the leading two terms.
The parameters $a_{0}$ and $r_{0}$ are determined such that the relation Eq.
(14) is optimized simultaneously for all the considered $D\bar{D}^{*}$ energy
levels. Our procedure follows the determinant residual method described in
Morningstar  _et al._ (2017), while additional details on implementation can
be found in appendix C of Ref. Prelovsek  _et al._ (2021).

We present two cases for each $J^{PC}=$ $1^{+-}$ and $1^{++}$: one where the
full correlation matrices are employed and the other where the meson-meson
interpolators with $\rho$ and/or $a_{0}$ mesons are omitted from the basis.
Therefore, four different fits are performed. The resulting scattering lengths
$a_{0}$ and the effective ranges $r_{0}$ are positive, as presented in Table 3
and Figs. 6, 7. Four and five $D\bar{D}^{*}$ energy levels from different
irreps are incorporated in the fits for $1^{+-}$ and $1^{++}$, respectively
(see the legends in the figures). We have verified that including additional
energy levels has a negligible effect on the resulting fit parameters.

Specific energy levels exhibit a negative shift, indicating attraction. This
results in a positive $k\cot{\delta_{0}}$ that increases with $k^{2}$.
However, some levels are consistent with noninteracting energies, resulting in
significant errors in the values of $k\cot{\delta_{0}}$ and consequently large
errors in the fit parameters. Slightly below the threshold, the central values
of these fits intersect with the dashed cyan line representing $ik=+|k|$ (see
the magenta octagons). This gives us a pole in the scattering amplitude (12)
at a real energy below the $D\bar{D}^{*}$ threshold
($E_{\textrm{cm}}^{\textrm{p}}<m_{D}+m_{D^{*}}$). Specifically, the poles
correspond to a virtual state, located at the scattering momentum
$k_{\textrm{p}}$, where $k_{\textrm{p}}=-i|k_{\textrm{p}}|$.101010The
scattering amplitude below the virtual state pole is not constrained by our
eigenenergies related to this channel, and we do not attribute physical
significance to $t_{\ell=0}$ at those energies.

Table 3: Table with parameters of four different fits to the $D\bar{D}^{*}$ scattering amplitude in the one-channel approximation; the scattering length, effective range and pole energy of the virtual state $\Delta m_{\textrm{V}}=E_{\textrm{cm}}^{\textrm{p}}-m_{D}-m_{D^{*}}$. Two fits for each $J^{PC}=1^{+-}$ and $1^{++}$ are performed, one from the analysis with full correlation matrices and the other one with certain interpolators excluded, as discussed in the text. $J^{PC}$ | interpolators |  | $1/a_{0}\,[\textrm{fm}^{-1}]$
---
$r_{0}\,[\textrm{fm}]$ | $\chi^{2}/N_{\textrm{dof}}$ | $\Delta m_{\textrm{V}}\,[\textrm{MeV}]$
$1^{+-}$ | all | $0.46^{+1.16}_{-0.45}$ | $0.96^{+0.43}_{-0.73}$ | $0.13$ | $-3.0^{+3.0}_{-31.1}$
$\eta_{c}\rho\ \textrm{excl.}$ | $0.54^{+1.07}_{-0.44}$ | $2.23^{+0.95}_{-1.08}$ | $0.24$ | $-2.8^{+2.6}_{-17.1}$
$1^{++}$ | all | $0.62^{+1.30}_{-0.51}$ | $1.78^{+0.25}_{-2.44}$ | $0.18$ |  $-3.8^{+3.6}_{}$111111Uncertainty is so large that it is unbounded from below.
$J/\psi\rho,\,\eta_{c}a_{0}\ \textrm{excl.}$ | $0.96^{+1.42}_{-0.91}$ | $2.19^{+0.36}_{-1.00}$ | $0.15$ | $-6.7^{+6.7}_{-19.5}$

A pole below the threshold can enhance the scattering rate above the
threshold, where the enhancement is more significant if the pole is closer to
the threshold. The effect of our virtual poles on the s-wave $D\bar{D}^{*}$
scattering rate (which is proportional to $k\left|t_{0}\right|^{2}$) is shown
in Fig. 8. Due to the large uncertainties in the scattering amplitude, the
significance of the enhancement is similarly uncertain.

A virtual state pole slightly below the threshold was also found in a recent
EFT study Zhang  _et al._ (2024), which argues why the molecular state with
$1(1^{++})$ does not produce a sizable experimental signal and has not been
observed so far. $Z_{c}$ has also been found as a virtual pole in a
phenomenological study of $D\bar{D}^{*}$ scattering Guo  _et al._ (2013).

The standard Lüscher scattering formalism breaks down when applied to energies
below the left-hand-cut (lhc) arising from the single-pion exchange. It
predicts a real scattering amplitude in cases where it should be complex. We
have not considered the effects stemming from the lhc, which is—in our
scenario with an unphysically high pion mass—at
$\left(ak_{\textrm{lhc}}\right)^{2}<-0.003$. This renders the extracted one-
channel phase shifts below the lhc unreliable and an effective range expansion
of the scattering amplitude is insufficient. Addressing these issues would
require a modification of the standard Lüscher method Dawid  _et al._ (2023);
Du  _et al._ (2023); Raposo and Hansen (2023); Meng  _et al._ (2024); Hansen
_et al._ (2024).

The results obtained under the strong assumptions within this section should
be taken as a first step to unveil the hadronic features hidden behind the
densely populated finite-volume eigenenergies we extract. These results
qualitatively suggest only a feeble interaction between $D$ and $\bar{D}^{*}$,
which is not strong enough to provide a real bound state, while it is not
possible to make further conclusions at this stage. An obvious next step is to
consider extracting the energy dependence of the amplitude in the coupled
channel system, including at least one of the channels in the calculation that
is ignored in the above analysis. We discuss such an approach in the next
section.

##  VII Reconciling $Z_{c}(3900)$ data from lattice studies and experiment?

An innovative aspect of the recent study presented in Ref. Yan  _et al._
(2024) involves conducting a comprehensive analysis that combines both
experimental $J/\psi\pi$, $D\bar{D}^{*}$ event distributions
$\textrm{d}N/\textrm{d}E_{\textrm{cm}}$ Ablikim  _et al._ (2017a, 2015b) and
the lattice finite-volume energy levels determined in Refs. Cheung  _et al._
(2017); Chen  _et al._ (2019). In this section, we utilize the same approach
as in Yan  _et al._ (2024), including our main lattice results in the analysis
along with the results of the other works. Compared to Cheung  _et al._
(2017); Chen  _et al._ (2019), our lattice results provide additional
information, in particular, valuable energy levels from the system in the
moving frame and higher energies in the rest frame. This section only
considers the $J^{PC}=1^{+-}$ channel and addresses the
$J/\psi\pi$-$D\bar{D}^{*}$ coupled system within a covariant effective field
theory framework. The aim is to constrain the unknown parameters in the
$J/\psi\pi$ and $D\bar{D}^{*}$ scattering amplitudes, thereby providing more
definitive insights into the properties of $Z_{c}(3900)$. For a detailed
derivation of the fitted experimental event distributions
$\textrm{d}N/\textrm{d}E_{\textrm{cm}}$, the procedure to relate scattering
amplitudes and finite-volume energies, the parameters used and other details,
we refer the reader to Ref. Yan  _et al._ (2024). Here, we only briefly
summarize the approach.

###  VII.1 The underlying EFT and the procedure behind the global fits

We employ the effective Lagrangians within a relativistic framework,
specifically emphasizing the energy region near the $D\bar{D}^{*}$ threshold.
As argued later, no pion exchange terms are considered, so interactions are
approximated using only local contact four-meson terms, guided by the
symmetries inherent in the theory:

  * •

The diagonal coupling $D\bar{D}^{*}\leftrightarrow D\bar{D}^{*}$ is obtained
from the most general Lagrangian consistent with the heavy quark spin symmetry
incorporated via
$H=\left[D^{*}_{\mu}\gamma^{\mu}-D\gamma_{5}\right](1+\not{v})/2$ AlFiky  _et
al._ (2006); Mehen and Powell (2011)

$\displaystyle\mathcal{L}_{D\bar{D}^{*}D\bar{D}^{*}}=$ $\displaystyle-
C_{0a}\left(D^{\dagger}D+D_{\mu}^{*\dagger}D^{*\mu}\right)\left(\bar{D}\bar{D}^{\dagger}+\bar{D}_{\mu}^{*}\bar{D}^{*\mu\dagger}\right)$
(15)
$\displaystyle+C_{0b}\left(D^{\dagger}D_{\mu}^{*}+D_{\mu}^{*\dagger}D\right)\left(\bar{D}^{*\mu}\bar{D}^{\dagger}+\bar{D}\bar{D}^{*\mu\dagger}\right)\>.$

Only terms at the leading-order $\mathcal{O}(p^{0})$ in the chiral expansion
are retained, while next-to-leading-order terms $\mathcal{O}(p^{2})$ are
neglected. The latter have a minor contribution since the heavy meson momenta
are small. For the $1^{+-}$ channel, only one combination of $C_{0a}$ and
$C_{0b}$ contributes Guo  _et al._ (2013); AlFiky  _et al._ (2006); Pavon
Valderrama (2012) and we refer to it as $\hat{\lambda}_{1}$.

Pion exchange interaction is not included explicitly. One reason is that the
short-range pion exchange can be absorbed into the redefinition of contact
interactions employed. Also, the coupling of pions to a pair of heavy mesons
involves a derivative, and therefore, this contribution is significantly
momentum suppressed near the threshold (at least at tree level). This
suppression is lifted at higher energies, which may be one of the reasons that
this EFT and lattice data show some tensions when including the data points
above the $D\bar{D}^{*}$ threshold, as will be discussed.

  * •

Regarding the contact interactions between $J/\psi\pi$ and $D\bar{D}^{*}$, the
general covariant operators featuring the fewest derivatives and being
invariant under $C$, $P$, chiral, and isospin symmetry transformations are
described by

$\displaystyle\mathcal{L}_{D\bar{D}^{*}J/\psi\pi}$
$\displaystyle=\hat{\lambda}_{2}\psi_{\mu}\left(\nabla^{\mu}D^{\dagger}u_{\nu}\bar{D}^{*\nu\dagger}+\bar{D}^{*\nu}u_{\nu}\nabla^{\mu}D\right)$
(16)
$\displaystyle+\hat{\lambda}_{3}\psi_{\mu}\left(\nabla^{\nu}D^{\dagger}u^{\mu}\bar{D}_{\nu}^{*\dagger}+\bar{D}_{\nu}^{*}u^{\mu}\nabla^{\nu}D\right)$
$\displaystyle+\hat{\lambda}_{4}\psi_{\mu}\left(\nabla^{\nu}D^{\dagger}u_{\nu}\bar{D}^{*\mu\dagger}+\bar{D}^{*\mu}u_{\nu}\nabla^{\nu}D\right)$
$\displaystyle+\hat{\lambda}_{5}\psi_{\mu}\left(D^{\dagger}\nabla^{\mu}u^{\nu}\bar{D}_{\nu}^{*\dagger}+\bar{D}_{\nu}^{*}\nabla^{\mu}u^{\nu}D\right)\>.$

Here, $\psi_{\mu}$ is the $J/\psi$ field, whereas $u_{\mu}$ refers to the
chiral pseudoscalar meson fields as defined in Refs. Gong  _et al._ (2016);
Yan  _et al._ (2024). This Lagrangian was used in Ref. Gong  _et al._ (2016)
and initially in a nonrelativistic form in Mehen and Powell (2011).

  * •

The $\mathcal{L}_{J/\psi\pi J/\psi\pi}$ is neglected since the interaction
between a pion and a heavy quarkonium is suppressed. This agrees with the
small scattering length of the $J/\psi\pi$ interaction (see Sec. IX and Refs.
Yokokawa  _et al._ (2006); Liu  _et al._ (2008); Liu (2009); Liu  _et al._
(2013b)) and is discussed also in Hanhart  _et al._ (2015); Guo  _et al._
(2016).

The Lagrangians (15) and (16) render the transition amplitudes for
$\bar{D}^{*}(a)D(b)\to\bar{D}^{*}(c)D(d)$ and
$J/\psi(a)\pi(b)\to\bar{D}^{*}(c)D(d)$,

$V_{\bar{D}^{*}D\rightarrow\bar{D}^{*}D}=\lambda_{1}\varepsilon_{c}^{\dagger}\cdot\varepsilon_{a}$
(17)

and

$\displaystyle V_{J/\psi\pi\rightarrow\bar{D}^{*}D}=\frac{1}{F_{\pi}}\Big{(}$
$\displaystyle\lambda_{2}\left(\varepsilon_{a}\cdot
p_{d}\right)\left(\varepsilon_{c}^{\dagger}\cdot p_{b}\right)$ (18)
$\displaystyle+$ $\displaystyle\lambda_{3}\left(\varepsilon_{a}\cdot
p_{b}\right)\left(\varepsilon_{c}^{\dagger}\cdot p_{d}\right)$
$\displaystyle+$
$\displaystyle\lambda_{4}\left(\varepsilon_{c}^{\dagger}\cdot\varepsilon_{a}\right)\left(p_{b}\cdot
p_{d}\right)$ $\displaystyle+$
$\displaystyle\lambda_{5}\left(\varepsilon_{a}\cdot
p_{b}\right)\left(\varepsilon_{c}^{\dagger}\cdot p_{b}\right)\Big{)}\>,$

where parameters $\lambda_{j}$ are proportional to $\hat{\lambda}_{j}$ as
detailed in Yan  _et al._ (2024). 121212Notice there is a typo for Eqs. (6)
and (7) in Ref. Yan  _et al._ (2024): the factor of $\sqrt{2}$ should be one.

It is then convenient to proceed with partial-wave amplitudes. In scattering
processes with spinful particles, both the $\ell S$ and helicity bases can be
employed for the partial-wave projections. While these two methods are
generally equivalent, the $\ell S$ basis is better suited for our study. The
reason is that in the molecular picture of $Z_{c}(3900)$, the s-wave
interaction of $D\bar{D}^{*}$ should dominate, with the d-wave part
significantly suppressed, at least when focusing on the energy region near the
$D\bar{D}^{*}$ threshold. We adopt the approach from Ref. Gülmez  _et al._
(2017) to perform partial-wave projections in a covariant manner. This
naturally introduces specific energy-dependent terms to the amplitudes from
polarization vectors, which are higher-order in the nonrelativistic expansion,
without adding extra free parameters. Refs. Albaladejo  _et al._ (2016); Du
_et al._ (2022) emphasize the crucial role of energy dependence in the
interaction kernel ($V(s)$) for generating resonance poles of $Z_{c}(3900)$ in
the complex energy plane near the physical region. Therefore, the covariant
scattering amplitudes approach should enable us to gain important insights
into the properties of $Z_{c}(3900)$.

The $J/\psi\pi$ and $D\bar{D}^{*}$ channels are labeled as channels 1 and 2,
respectively. The s-wave amplitudes $V_{11}(s)$131313In the last bullet point
in this section, we have already mentioned that the perturbative $J/\psi\pi\to
J/\psi\pi$ transition amplitude is insignificantly small. Consequently, the
corresponding matrix element vanishes, $V_{11}=0$., $V_{12}(s)$141414The
$\hat{\lambda}_{3}$ term in Eq. (16) has no impact on the s-wave amplitude.,
$V_{22}(s)$ build up a symmetric matrix spanned in the scattering-channel
space,

$V(s)=\left[\begin{array}[]{ll}V_{11}(s)&V_{12}(s)\\\
V_{12}(s)&V_{22}(s)\end{array}\right]\>,$ (19)

where $\sqrt{s}=E_{\textrm{cm}}$ and the explicit expressions of the partial-
wave amplitudes $V_{ij}$ can be found in Ref. Yan  _et al._ (2024). $V(s)$ is
further used to derive the on-shell unitary partial-wave two-body scattering
amplitude151515This scattering amplitude and the one defined in Sec. VI are
related, with $T=-16\pi t$.

$T(s)=[1-V(s)\cdot G(s)]^{-1}\cdot V(s)\>,$ (20)

where $G(s)$ is the well-known diagonal loop function matrix that also
contains the two free parameters $a_{\textrm{SC,}i}$, $i=1,2$, called
subtraction constants. This equation is diagrammatically shown in Fig. 9.

Figure 9: Diagrammatic representation of Eq. (20), which can be rewritten as
$T=V+VGT$, where $G$ denotes the loop.

The experimental $J/\psi\pi$ and $D\bar{D}^{*}$ invariant-mass distributions

$\frac{dN_{i}}{d\sqrt{s}}=A_{i}(P_{i}(s))+B_{i}(s),\quad i=1,2\>$ (21)

are projected from the three-body decays $Y\to J/\psi\pi\pi$ and $Y\to
D\bar{D}^{*}\pi$, respectively. $A_{i}$ and $B_{i}$ represent the signal and
background part, respectively. The former is a function of the two-body
production amplitudes (see Fig. 10)

$\footnotesize{\mathcal{P}(s)=\left[\begin{array}[]{l}P_{1}(s)\\\
P_{2}(s)\end{array}\right]=[1-V(s)\cdot
G(s)]^{-1}\cdot\alpha\>,\quad\alpha=\left[\begin{array}[]{c}\alpha_{1}\\\
\alpha_{2}\end{array}\right]\approx\left[\begin{array}[]{c}0\\\
\alpha_{2}\end{array}\right]\,}$ (22)

where $\alpha$ are constant production vertices. Given that the cross section
of $e^{+}e^{-}\to\pi^{+}D^{0}D^{*-}$ between $4.2\text{\,}$ and
$4.3\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ is much larger
(200–$300\text{\,}\mathrm{p}\mathrm{b}$ Ablikim  _et al._ (2019)) than that of
$e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ (50–$80\text{\,}\mathrm{p}\mathrm{b}$
Ablikim  _et al._ (2017b)), we set the direct production-vertex constant
$\alpha_{1}$ to zero. This is also consistent with the choice of setting
$V_{11}=0$ in Eq.(19).

Figure 10: Diagrammatic representation of the production amplitude in Eq.
(22), which can be rewritten as $\mathcal{P}=\alpha+TG\alpha$, where $G$
represents the loop. The negligible $\alpha_{1}$ is represented with the faded
text of the direct product channel $J/\psi\pi$.

We will perform a global fit to the experimental and lattice data. Once we
establish the values of the unknown parameters in the production amplitudes
(22), the unitarized scattering amplitudes (20) become unambiguously
determined. Consequently, we can extract information about the $Z_{c}(3900)$
resonance, including its pole positions and coupling strengths in the relevant
channels.

To fit also the lattice data, a finite box of length $L$ with periodic
boundary conditions and inertial frame with total momentum
$|{\hbox{\boldmath{$P$}}}|=0$ or $2\pi/L$ is considered. The potential $V(s)$
and the finite volume energies $E_{\textrm{cm}}=\sqrt{s}$ observed on the
lattice are related via the Lüscher’s quantization condition, which was
conveniently expressed in Refs. Döring  _et al._ (2011); Doring  _et al._
(2012); Guo  _et al._ (2019); Gockeler  _et al._ (2012) as

$\det\left[1-V(s)\cdot\left(G(s)+\Delta G(s)\right)\right]=0\>~{}.$ (23)

Here, $\Delta G(s)$ is the finite volume correction of the loop function. We
remark that correcting our lattice energies with (9) and (10) enables our
lattice energies to be jointly fitted with energies from Refs. Cheung  _et
al._ (2017); Chen  _et al._ (2019), where a relativistic dispersion relation
is employed. This is a notable improvement over the approach in Ref. Prelovsek
_et al._ (2015). Comparing with Ref. Yan  _et al._ (2024), where only lattice
data in the rest frame were included, the present work also considers the EFT
fits to the lattice energy levels in the moving frame. We refer to Ref. Guo
_et al._ (2019), especially Eqs. (31–35) of this reference, for the details
about calculating the moving-frame spectra in the unitarized EFT approach.

###  VII.2 Global fits to the experimental and lattice data

The unknown parameters are determined by simultaneously fitting experimental
and lattice data, as in Ref. Yan  _et al._ (2024). Two types of relevant
experimental data are included: first, the $J/\psi\pi^{\pm}$ event
distributions from the $e^{+}e^{-}\to J/\psi\pi^{+}\pi^{-}$ process at center-
of-momentum energies $4.23\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and
$4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ Ablikim  _et al._ (2017a), and
second, the $D^{0}D^{*-}$ and $D^{-}D^{*0}$ event distributions from the
$e^{+}e^{-}\to\pi^{\pm}(D\bar{D}^{*})^{\mp}$ processes at the same
$e^{+}e^{-}$ energies Ablikim  _et al._ (2015b). For lattice data, we take the
finite-volume $J/\psi\pi$ and $D\bar{D}^{*}$ energies from the left-hand sides
of Figs. 1, 2, and data from Refs. Cheung  _et al._ (2017); Chen  _et al._
(2019). The fits employ the masses of the scattered hadrons that feature in
the lattice simulations (see Appendix A and Ref. Yan  _et al._ (2024)).

We assume that the fitted parameters do not depend on masses of scattering
hadrons and then provide the scattering amplitudes and poles for the case of
physical masses. Our parameters are the dimensionless coupling constants
$\lambda_{1}$ and
$\tilde{\lambda}_{j=2,4,5}\equiv\lambda_{j=2,4,5}\,m_{D}^{\rm phys}$
161616There is a typo for Eq. (24) in Ref. Yan  _et al._ (2024): it should be
$\tilde{\lambda}_{j=2,3,4,5}=\lambda_{j=2,3,4,5}\,m_{D}^{\rm phys}$.
($m_{D}^{\rm phys}$ being the physical $D$-meson mass), two background
parameters $b_{1}$ from the $Y\to J/\psi\pi\pi$ processes at both
$4.23\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and
$4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ center-of-momentum energies, the
four production-vertex constants $\alpha_{2}$ for both $e^{+}e^{-}\to
J/\psi\pi^{+}\pi^{-}$ and $e^{+}e^{-}\to\pi^{\pm}(D\bar{D}^{*})^{\mp}$
processes at both energies and the subtraction constants $a_{\textrm{SC,}1}$,
$a_{\textrm{SC,}2}$. The background events of the $D\bar{D}^{*}$ event
distribution are subtracted from the analysis of Ref. Ablikim  _et al._
(2015b). Consequently, our fitting approach excludes the background term
$B_{2}(s)$ from (21) but retains $B_{1}(s)$, which includes the parameter
$b_{1}$. The subtraction constants, introduced as part of the unitarization
process, are left to vary freely.

Four distinct fits—all including experimental and lattice data—are conducted,
and the reason for this will become evident later. The fits incorporate the
lattice levels relevant for the coupled scattering of $D\bar{D}^{*}$ and
$J/\psi\pi$ in s-wave. The common feature of all available lattice studies is
that all energies are close to noninteracting energies. The distinction
between the fits lies in the handling of four of our energy levels with
$E_{\textrm{cm}}>$4050\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$:

  * •

Fit1 is the conventional and preferred fit that includes all the data: the
experimental data, all 18 relevant lattice energy levels from our
simulation171717Note that no fit incorporates states with significant overlap
with $D(1)\bar{D}^{*}(2)$ and $D(2)\bar{D}^{*}(1)$ operators. Their energies
lie around $200\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ above the lattice
$D\bar{D}^{*}$ threshold, which is the edge of the energy region accessed by
the experiment. This and the fact that such high energy levels have not been
reliably determined represent the arguments for not including the
aforementioned energies. (red and green circles in Fig. 12) and levels from
simulations Cheung  _et al._ (2017); Chen  _et al._ (2019)181818Fits here and
in Yan  _et al._ (2024) incorporate the lowest six levels from Chen  _et al._
(2019)..

  * •

Fit2: same as Fit1, but the four of our levels at highest energies (green in
Fig. 12) are not included.

  * •

Fit3: same as Fit1, but our four highest levels are included with errors that
are artificially reduced by a factor of four (variations of this fit render
similar conclusions as discussed in Appendix F).

  * •

Fit4: the original joint fit from Ref. Yan  _et al._ (2024), which does not
incorporate our lattice levels.

So, Fit2 excludes four of our highest energy levels, while Fit3 gives these
levels more weight in the fit.

Table 4: Fitted parameters and the corresponding $\chi^{2}/N_{\textrm{dof}}$. All four fits include data from experiment Ablikim  _et al._ (2017a, 2015b) and previous lattice calculations Cheung  _et al._ (2017); Chen  _et al._ (2019). Additionally, Fit2 includes 14 data points from our new lattice simulation below $4050\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$, while Fit1 and Fit3 include 18 of our data points below $4150\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$. For Fit3, we have artificially reduced the uncertainties of our four lattice levels between $4050\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ and $4150\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ by a factor of 4. The left/right values for the $b_{1}$ and $\alpha_{2}$ represent results for the data sets at center-of-momentum energies $4.23\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and $4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, respectively.  parameter |  Fit1 (preferred) |  Fit2 |  Fit3 |  Fit4 Yan  _et al._ (2024)
---|---|---|---|---
$a_{\textrm{SC},1}$ |  $-6.63^{+1.87}_{-3.69}$ |  $-4.44^{+0.43}_{-0.52}$ |  $-9.27^{+2.87}_{-4.23}$ |  $-4.02^{+0.32}_{-0.49}$
$a_{\textrm{SC},2}$ |  $-2.80^{+0.06}_{-0.17}$ |  $-2.81^{+0.03}_{-0.03}$ |  $-2.05^{+0.07}_{-0.11}$ |  $-2.80^{+0.02}_{-0.02}$
$\lambda_{1}$ |  $-102^{+36}_{-13}$ |  $-88^{+11}_{-12}$ |  $-63^{+71}_{-56}$ |  $-86^{+12}_{-11}$
$\tilde{\lambda}_{2}$ |  $219^{+151}_{-172}$ |  $999^{+115}_{-157}$ |  $-2208^{+602}_{-489}$ |  $1082^{+93}_{-116}$
$\tilde{\lambda}_{4}$ |  $36^{+9}_{-13}$ |  $-44^{+17}_{-11}$ |  $340^{+71}_{-89}$ |  $-41^{+14}_{-11}$
$\tilde{\lambda}_{5}$ |  $160^{+119}_{-133}$ |  $-734^{+185}_{-121}$ |  $3399^{+710}_{-883}$ |  $-751^{+137}_{-111}$
$b_{1}\left[10^{-4}\times\rm{MeV^{-3}}\right]$ |  $9.0^{+0.2}_{-0.2}/4.9^{+0.1}_{-0.1}$ |  $8.8^{+0.2}_{-0.2}/4.8^{+0.1}_{-0.1}$ |  $8.8^{+0.2}_{-0.2}/4.8^{+0.1}_{-0.1}$ |  $8.7^{+0.2}_{-0.2}/4.8^{+0.1}_{-0.2}$
$\left|\alpha^{J/\psi\pi}_{2}\right|^{2}$ |  $30.8^{+6.3}_{-10.0}/14.1^{+2.9}_{-4.7}$ |  $20.2^{+3.4}_{-2.5}/9.3^{+1.8}_{-1.2}$ |  $316.4^{+127.2}_{-117.7}/147.3^{+56.6}_{-55.5}$ |  $19.4^{+3.4}_{-2.1}/9.0^{+1.6}_{-1.1}$
$\left|\alpha^{D\bar{D}^{*}}_{2}\right|^{2}$ |  $1.3^{+0.5}_{-0.4}/0.6^{+0.2}_{-0.2}$ |  $2.2^{+0.5}_{-0.4}/1.1^{+0.2}_{-0.2}$ |  $6.7^{+1.4}_{-1.6}/2.8^{+0.6}_{-0.7}$ |  $2.7^{+0.6}_{-0.6}/1.3^{+0.3}_{-0.3}$
$\chi^{2}/N_{\textrm{dof}}$ |  $\frac{484.1}{368-12}=1.36$ |  $\frac{412.8}{364-12}=1.17$ |  $\frac{645.8}{368-12}$=1.81 |  $\frac{398.7}{350-12}=1.18$

The fitted parameters of all four fits are represented in Table 4. The
parameters of Fit2 agree very well with Fit4. On the other hand, certain
parameters of Fit1 and Fit3 are not consistent with Fit2 and Fit4. This is
more pronounced for Fit3. The $\chi^{2}/N_{\mathrm{dof}}$ for Fit1 and Fit3 is
larger than for the other two. However, the resulting
$\chi^{2}/N_{\mathrm{dof}}$ for Fit1 is still relatively small.

Predictions from the fits are compared to the experimental invariant mass
distributions for $D\bar{D}^{*}$ and $J/\psi\pi$ in Fig. 11. They show good
agreement for both $e^{+}e^{-}$ energies. In particular, peaks related to the
$Z_{c}$ feature for all fits. Only Fit3 somewhat underestimates the height of
the peak in the $D\bar{D}^{*}$ distribution.

Figure 11:  The $J/\psi\pi$ (top) and $D\bar{D}^{*}$ (bottom) event
distributions from experiment (symbols) and from the EFT approach (lines). The
top pane shows $J/\psi\pi^{\pm}$ invariant mass distributions of the process
$e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ at
$E_{e^{+}e^{-}}=$4.23\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ and
$4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ from BESIII Ablikim  _et al._
(2017a). The bottom pane shows the $D^{0}D^{*-}$ and $D^{-}D^{*0}$ event
distributions of $e^{+}e^{-}\to\pi^{\pm}(D\bar{D}^{*})^{\mp}$ from BESIII
Ablikim  _et al._ (2015b) at the same $E_{e^{+}e^{-}}$. Background events from
the experimental analysis were subtracted. Blue, orange and green lines
represent Fit1, Fit2 and Fit3, respectively. The shaded areas are the
uncertainties of Fit1.

Let us confront predictions of the fits with lattice energies from the present
and previous simulations: the HadSpec collaboration provides energy levels up
to the $D\bar{D}^{*}$ threshold Cheung  _et al._ (2017), while CLQCD Chen  _et
al._ (2019) and the present study also provide several levels above the
threshold. The common feature is that all energies are close to the
noninteracting energies. Our lattice energies are compared with the preferred
Fit1 in Fig. 12. All fits are compared to our and other lattice data in Figs.
13 and 14. Note that one should not directly compare the energies from
different simulations since they employ different quark masses and volumes.

The predicted energies of the preferred conventional Fit1 manifest some
disagreement with our lattice levels above
$E_{\textrm{cm}}=$4050\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ and we
elaborate on possible effects that cause this in section VII.6. This
discrepancy is also the motivation for performing several fits. The high-lying
energy levels are not included in Fit2, which achieves a smaller
$\chi^{2}/N_{\mathrm{dof}}$ than Fit1. On the other hand, Fit3 nearly
reproduces our high-lying energy levels, although they lie close to the
noninteracting energies. In summary, the three different fits produce
different results since they treat the four decisive lattice data points
differently. Fit3 artificially emphasizes these data points, Fit2 excludes
them, while the conventional Fit1 includes them with correct weights. This is
the reason why we have chosen Fit1 as the preferred fit.

Figure 12:  Energy levels $E_{\textrm{cm}}$ (red and green squares) from our
lattice calculations together with the prediction of the preferred Fit1 (blue
solid lines) as a function of $L$. This fit incorporates data points from
experimental measurements Ablikim  _et al._ (2017a, 2015b), from lattice
studies Cheung  _et al._ (2017); Chen  _et al._ (2019) and 18 of our lattice
data points. Total momentum is $\left|{\hbox{\boldmath{$P$}}}\right|=0$ on the
left and $\left|{\hbox{\boldmath{$P$}}}\right|=1\cdot 2\pi/L$ on the right
panel. The dark red dashed lines denote the free energy levels. The gray
shadowed areas represent the error bands from Fit1. Only the finite-volume
energy levels that are dominantly coupled to s-wave $J/\psi\pi$ and
$D\bar{D}^{*}$ interpolators are shown and used in the fit.

Figure 13:  Lattice energy levels $E_{\textrm{cm}}$ (red, green and brown
squares) for zero total momentum $\left|{\hbox{\boldmath{$P$}}}\right|=0$ from
our main lattice results (left), results from HSC (right) Cheung  _et al._
(2017) and CLQCD (in the middle) Chen  _et al._ (2019). They are shown
together with two different types of performed fits, Fit2 (orange dotted-
dashed lines) and the preferred Fit1 (blue solid line). Both fits incorporate
data points from experimental measurements Ablikim  _et al._ (2017a, 2015b)
and from lattice studies Cheung  _et al._ (2017); Chen  _et al._ (2019). Fit2
includes 14 of our lattice data points (red squares), while Fit1 includes 18
(red and green squares). The highest brown CLQCD energy point was not used in
any of the fits. The dark red dashed lines denote the free energy levels. The
gray shadowed areas represent the error bands from Fit1. We omit the error
bands from Fit2 to avoid overloading the plots. Only the finite-volume energy
levels that are dominantly coupled to s-wave $J/\psi\pi$ and $D\bar{D}^{*}$
interpolators are shown and used in the fit.

Figure 14:  As Fig. 13, but Fit3 (green dotted-dashed lines) is shown instead
of Fit2. The uncertainties of four high-lying levels (shown by green squares)
are artificially reduced by a factor of four in Fit3.

The coupling constants $\lambda_{1}$ and $\tilde{\lambda}_{j=2,4,5}$ as well
as two subtraction constants characterize the nonperturbative interactions of
$J/\psi\pi$ and $D\bar{D}^{*}$. These parameters are crucial in the coupled-
channel scattering amplitudes, impacting event distributions across all energy
points and influencing the lattice finite-volume spectra. These six parameters
simultaneously influence both the experimental and lattice data under
consideration. In contrast, the background parameters $b_{1}$ and the
production vertex parameters $\alpha_{2}$, which serve as normalization
constants, are anticipated to vary with the different BESIII datasets at
$4.23\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and
$4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$.

###  VII.3 Insights into the $Z_{c}(3900)$ resonance poles

Kinematical singularities, e.g., two-body cusps and triangle and box
singularities, offer a potential explanation for specific experimentally
observed structures Guo  _et al._ (2020). These singularities exhibit high
sensitivity to the specific kinematics of the processes. In contrast,
resonance poles in the complex energy plane possess a universality that
extends across all production amplitudes involving the same particles. The
agreement of the employed model with experimental invariant mass distributions
at both $e^{+}e^{-}$ energies suggests that the experimental $Z_{c}$ peak does
not decisively arise from the triangular singularity but rather from the
pole(s) Wang  _et al._ (2013); Albaladejo  _et al._ (2016); Du  _et al._
(2022); Gong  _et al._ (2018); Guo  _et al._ (2015). Therefore, we next search
for relevant poles in the system under investigation.

This section focuses only on the poles around the $D\bar{D}^{*}$ threshold. In
the vicinity of a pole singularity at $s=s_{\textrm{R}}$, the elements of the
scattering amplitude follow a behavior described by:

$T_{ij}(s)=-\frac{\gamma_{i}\gamma_{j}}{s-s_{\textrm{R}}}\>,$ (24)

where indexes $i$ and $j$ stand for the (in this case, two) coupled channels.
The real and imaginary parts of the pole position are conventionally
identified as the mass and half-width of a resonance, represented as
$\sqrt{s_{\textrm{R}}}=m_{\textrm{R}}\pm\frac{i}{2}\Gamma_{\textrm{R}}$. The
complex-valued residues $\gamma_{i}$ offer insights into the strength of the
resonance’s coupling to scattering channels.

The complex energy plane exhibits branch cuts associated with the opening of
new channels. For single-channel scattering, the physical and unphysical
Riemann sheets (RS) correspond to $\textrm{Im}(k)>0$ and $\textrm{Im}(k)<0$,
respectively. Here $k$ is the magnitude of the three-momentum in the cm-frame.
In our two-channel case ($J/\psi\pi$, $D\bar{D}^{*}$), the following labeling
of sheets is adopted:

RS I:
$\displaystyle\qquad\textrm{Im}(k_{J/\psi\pi})>0\quad\textrm{Im}(k_{D\bar{D}^{*}})>0\>,$
(25) RS II:
$\displaystyle\qquad\textrm{Im}(k_{J/\psi\pi})<0\quad\textrm{Im}(k_{D\bar{D}^{*}})>0\>,$
RS III:
$\displaystyle\qquad\textrm{Im}(k_{J/\psi\pi})<0\quad\textrm{Im}(k_{D\bar{D}^{*}})<0\>,$
RS IV:
$\displaystyle\qquad\textrm{Im}(k_{J/\psi\pi})>0\quad\textrm{Im}(k_{D\bar{D}^{*}})<0\>.$

The physical measurements are performed along the real energy axes above
threshold on the physical Riemann sheet I.

Two poles relatively close to $D\bar{D}^{*}$ threshold and to the physical
region are found for all fits. They are located on RS III and RS IV for the
central values of the preferred Fit1. Both poles seem to have comparable
influence according to their closest connections to the physical axis. This
can be seen in Figs. 15 and 16 containing pole singularities of the preferred
Fit1 in the complex energy and momentum planes. The resonance poles’ positions
from our fits and the corresponding residues are presented in Table 5. We
discuss only the poles near the real axis and the $D\bar{D}^{*}$ threshold.
Every bootstrap sample contains a resonance pole on RS III and another either
on RS II or RS IV, as in Yan  _et al._ (2024).

Table 5: Relevant poles of the $Z_{c}(3900)$ and their residues in different fits. These values have been obtained using the physical masses of the scattering particles. $\gamma_{1}$ and $\gamma_{2}$ are residues of the $J/\psi\pi$ and $D\bar{D}^{*}$ channels, respectively. All fits have two near-threshold poles. Since the error determination procedure involves bootstrap sampling, different samples can give poles on different RSs. A sample contains either poles on RS III and II or poles on RS III and IV (see Fig. 17).  RS |  $\quad m_{\textrm{R}}\,$[MeV] |  $\Gamma_{\textrm{R}}/2\,$[MeV] |  $\left|\gamma_{1}\right|\,$[GeV] |  $\left|\gamma_{2}\right|\,$[GeV]
---|---|---|---|---
Fit1 |  |  |  |
III | $3878.4^{+4.2}_{-3.4}$ | $31.9^{+3.3}_{-5.8}$ | $3.2^{+0.5}_{-0.5}$ | $10.6^{+0.5}_{-1.6}$
IV${}_{\mathrm{(central\ v.)}}$ | $3894.3^{+2.7}_{-3.2}$ | $20.9^{+3.2}_{-7.1}$ | $3.2^{+0.4}_{-0.6}$ | $10.1^{+0.4}_{-1.1}$
or |  |  |  |
II | $3900.8^{+4.1}_{-4.1}$ | $1.58^{+1.56}_{-1.56}$ | $4.7^{+0.6}_{-0.6}$ | $8.2^{+0.9}_{-0.9}$
Fit2 |  |  |  |
III | $3876.0^{+3.5}_{-4.6}$ | $32.0^{+2.5}_{-1.6}$ | $4.1^{+0.3}_{-0.3}$ | $9.1^{+1.0}_{-0.7}$
IV${}_{\mathrm{(central\ v.)}}$ | $3900.8^{+1.9}_{-1.9}$ | $7.7^{+5.4}_{-4.1}$ | $4.3^{+0.3}_{-0.4}$ | $8.9^{+0.6}_{-0.4}$
or |  |  |  |
II | $3902.2^{+1.2}_{-1.2}$ | $1.44^{+1.42}_{-1.42}$ | $4.8^{+0.3}_{-0.3}$ | $8.4^{+0.3}_{-0.3}$
Fit3 |  |  |  |
III | $3820.9^{+28.9}_{-52.6}$ | $41.5^{+23.0}_{-17.3}$ | $4.7^{+0.6}_{-0.8}$ | $16.2^{+4.5}_{-3.8}$
IV | $3858.6^{+11.9}_{-5.9}$ | $29.4^{+4.6}_{-11.8}$ | $4.1^{+0.3}_{-0.7}$ | $15.4^{+1.7}_{-2.6}$
Fit4 |  |  |  |
III | $3874.8^{+3.7}_{-4.0}$ | $32.7^{+1.6}_{-1.9}$ | $4.3^{+0.3}_{-0.3}$ | $8.7^{+0.8}_{-0.7}$
IV${}_{\mathrm{(central\ v.)}}$ | $3902.4^{+1.1}_{-2.3}$ | $3.3^{+6.3}_{-1.8}$ | $4.6^{+0.2}_{-0.4}$ | $8.6^{+0.6}_{-0.2}$
or |  |  |  |
II | $3902.7^{+1.3}_{-1.3}$ | $3.0^{+2.4}_{-2.4}$ | $4.9^{+0.2}_{-0.2}$ | $8.3^{+0.3}_{-0.3}$
Figure 15: Different paths in the complex energy plane that smoothly connect
poles on RS III and RS IV (from Fit1) with the region where the physical
$D\bar{D}^{*}$ scattering occurs. The green line represents the real energy
axis. The circle represents the physical $D\bar{D}^{*}$ threshold, while the
$J/\psi\pi$ threshold is not shown since it is too distant.

Figure 16: Locations of the resonance poles of the scattering amplitude based
on EFT for the physical masses of the scattered mesons. Results from the
preferred conventional Fit1 are shown. There are two poles, one on RS III and
another one on either RS II or RS IV. Left panel: the complex energy plane.
The green line represents the real energy axis. Right panel: the complex plane
of the $D\bar{D}^{*}$-momentum in the cm-frame (complex-$k_{D\bar{D}^{*}}$
plane) divided into four quadrants for four RSs. The subscripts $l/u$
represent the lower/upper complex energy half-plane. The cool and warm colored
contours represent constant $\textrm{Re}\left(E_{\textrm{cm}}\right)$ and
$-2\textrm{Im}\left(E_{\textrm{cm}}\right)$ in MeV, respectively. They aid
visualization of the proximity of poles to physical scattering, which occurs
along the positive imaginary axis (below the $D\bar{D}^{*}$) and the positive
real axis (above the $D\bar{D}^{*}$ threshold). Small circles represent the
spread of poles according to all bootstrap samples. The physical
$D\bar{D}^{*}$ threshold is represented by the circle (in the origin), while
the $J/\psi\pi$ threshold (on the imaginary axis) is not shown since it is too
distant.

Figure 17: Locations of the resonance poles of the scattering amplitude from
Fit1 (circles), Fit2 (triangles) and Fit3 (squares). For more details, see
Fig. 16.

Poles of the scattering amplitude from Fit1, Fit2 and Fit3 are presented in
Fig. 17. Let us now focus on those from Fit1 and Fit2. RS II/IV poles lie
close to each other, which suggests that they are qualitatively the same. On
the other hand, the pole on RS II/IV and one on RS III are not so close. Given
that the mass of the experimentally observed $Z_{c}(3900)$
($\approx$3887\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$) Workman  _et al._
(2022) is in between the real parts of both poles’ complex energies, it is
plausible that the $Z_{c}(3900)$ peak is a manifestation caused by both poles.

One can see that the poles from Fit3 are further away from the physical
region. Their position is still close enough to produce the peaks in the line
shapes (see Fig. 11), however, with somewhat smaller peaks in the
$D\bar{D}^{*}$ distributions (see green lines).

We have checked that all the mentioned poles disappear by removing the off-
diagonal matrix elements in the interaction kernel ($V_{12}$). The comparison
of the magnitude of the $D\bar{D}^{*}\to D\bar{D}^{*}$ scattering amplitude by
including/excluding the $J/\psi\pi\to D\bar{D}^{*}$ contribution is given in
Fig. 18. This indicates that the coupling between both channels is likely
important for the existence of the $Z_{c}$ peak (see also Subsec. VIII.1).

Figure 18: The magnitude of the diagonal $D\bar{D}^{*}\to D\bar{D}^{*}$
scattering amplitude, with and without $V_{12}$. This suggests that the
coupling between both channels is important for the existence of the near-
threshold peak and poles.

The pole singularities on the unphysical RSs in the complex energy plane have
noticeable effects on the physical amplitudes, which can be characterized by
the phase shifts ($\delta_{1}$, $\delta_{2}$) and the inelasticity ($\eta$) in
the coupled-channel scattering matrix191919We omit the labels $\ell=0$ for the
phase shifts denoting the s-wave $J/\psi\pi$ and $D\bar{D}^{*}$ scattering.

$S=\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]\>.$ (26)

The predicted phase shifts and inelasticity for Fit1 are presented in Fig. 19.
Within the defined uncertainties, two branches of phase shifts emerge above
the $D\bar{D}^{*}$ threshold. Our analysis affirms that the upper branch of
$\delta_{2}$ and the lower branch of $\delta_{1}$ correspond to parameter
samples with a pole on the RS IV. The other two branches correspond to
parameter samples with a pole on the RS II. Therefore, the rise of
$\delta_{2}$ through $90^{\circ}$ is likely caused by the pole on RS IV and
not the one on RS II. We can not rule out that the pole on RS III also
influences the $Z_{c}$ peak. Therefore, it is plausible that poles on RS IV
and III can both influence the physical peak.

Figure 19: Phase shifts of $J/\psi\pi\to J/\psi\pi$ (left), $D\bar{D}^{*}\to
D\bar{D}^{*}$ (center) and the inelasticity (right) as defined in (26). Shown
are the results for Fit1 (blue solid line), Fit2 (orange dashed line) and Fit3
(green dotted-dashed line) with the error band from Fit1.

The $Z_{c}(3900)$ (and its strange partner $Z_{cs}(3985)$) is addressed in Du
_et al._ (2022), yielding a RS III pole with a position comparable to our
results from Fits 1, 2 and 4, as shown in Table 5. A similar pole position on
RS III is reported in a study considering the $\pi\pi$ event distributions
alongside the $J/\psi\pi$ and $D\bar{D}^{*}$ ones Chen  _et al._
(2023).202020In fact, thanks to the communications of Meng-Lin Du and Yun-Hua
Chen, another pole on RS IV is present in an almost identical location to the
one on RS III in both mentioned studies Du  _et al._ (2022); Chen  _et al._
(2023). However, such a pole is not highlighted since it is considered to have
a minor influence on the physical region compared to the RS III pole. To our
knowledge, only one study finds two poles at different locations and being
close to the threshold Zhou and Xiao (2015). Further details on the pole
positions identified by other studies are provided in Subsec. VIII.1.

###  VII.4 Ratio of partial decay widths

Based on the unitarized decay amplitudes, we estimate the ratio of
$Z_{c}(3900)$ partial decay widths into $D\bar{D}^{*}$ and $J/\psi\pi$
channels, $R\equiv\Gamma_{Z_{c}\to D\bar{D}^{*}}/\Gamma_{Z_{c}\to J/\psi\pi}$.
This involves performing phase space integrals within the energy range of
$E_{\textrm{cm}}=(3900\pm 35)\,$MeV,212121We verified that the change of the
energy range changes the final values by much less than their $1\sigma$
uncertainties. following the same energy range and methodology proposed in
Ref. Albaladejo  _et al._ (2016). When evaluating the ratio of partial decay
widths for $Z_{c}(3900)$, we exclude background contributions and set the
$\alpha^{J/\psi\pi}_{2}$ and $\alpha^{D\bar{D}^{*}}_{2}$ factors—accounting
for the normalizations of experimental event distributions—to unity. Here are
the resulting ratios from various fits:

$\displaystyle R$ $\displaystyle=6.1^{+2.0}_{-2.3}~{}\quad({\rm Fit1})\>,$
(27) $\displaystyle R$ $\displaystyle=2.4^{+1.0}_{-0.7}~{}\quad({\rm
Fit2})\>,$ $\displaystyle R$ $\displaystyle=12.2^{+3.2}_{-3.0}~{}\quad({\rm
Fit3})\>,$ $\displaystyle R$ $\displaystyle=1.9^{+0.9}_{-0.6}~{}\quad({\rm
Fit4})\>.$

As for other quantities, the agreement between Fit2 and Fit4 also manifests
for $R$. Their values are somewhat lower than the experimentally determined
$R_{\textrm{exp}}=$6.2\text{\,}$\pm$1.1\text{\,}$\pm$2.7\text{\,}$$ Ablikim
_et al._ (2014a). An interesting finding is that Fit1, with a somewhat larger
$\chi^{2}$, is in closer agreement with experiment. On the other hand, the
result from Fit3 is even larger. However, one has to be careful with any
conclusions concerning this ratio since it was measured only once by BESIII
back in 2014 Ablikim  _et al._ (2014a).

###  VII.5 Compositeness

We can estimate the composition of $Z_{c}(3900)$ based on our results. The
partial two-hadron compositeness coefficient $X_{i}$ is calculated following
Ref. Guo and Oller (2016). It corresponds to the probability of having the
two-body component from channel $i$ in the considered resonance state. Since
we focus on the region around the $D\bar{D}^{*}$ threshold, we derive only the
corresponding $X_{2}$ for three fits

$\displaystyle X_{2}^{\textrm{III}}=0.65^{+0.05}_{-0.15},\quad
X_{2}^{\textrm{IV}}=0.59^{+0.04}_{-0.09}\quad~{}(\textrm{Fit1}),$ (28)
$\displaystyle X_{2}^{\textrm{III}}=0.48^{+0.10}_{-0.07},\quad
X_{2}^{\textrm{IV}}=0.46^{+0.05}_{-0.04}\quad~{}(\textrm{Fit2}),$
$\displaystyle X_{2}^{\textrm{III}}=0.43^{+0.09}_{-0.06},\quad
X_{2}^{\textrm{IV}}=0.42^{+0.05}_{-0.02}\quad~{}(\textrm{Fit4}).$

We do not consider compositeness coefficients for the Fit3 since the
corresponding poles are significantly below the threshold energy, while the
way to calculate $X$ in Ref. Guo and Oller (2016) is valid only for resonances
above the threshold. The resulting coefficients suggest that elementary and
$D\bar{D}^{*}$ molecular Fock components are almost equally important. Fit1
slightly favors the molecular component, which is also consistent with the $R$
values in (27) since the $R$ for Fit1 is higher than for Fit2 and Fit4.

###  VII.6 Interpretations

The locations of both poles suggest that they can equally contribute to the
physical peak. The ordinary fits (Fit1 and Fit2) of lattice and experimental
data yield a satisfactory replication of the invariant mass distributions. The
lattice eigenenergies lie close to noninteracting two-meson energies. The fits
nicely reproduce lattice energies below and at the $D\bar{D}^{*}$ threshold,
while reproduction above the threshold is not optimal. This does not
significantly worsen the fit (in terms of the $\chi^{2}$) since most of the
fitted points are experimental, not from the lattice. Still, we tried to put
more emphasis on lattice data and performed fits involving the four decisive
lattice energies with artificially reduced uncertainties (see Fit3 as the
example). As expected, the prediction gets closer to those lattice energy
points. One might expect that this fit would not reproduce the $Z_{c}(3900)$
peaks, but it does: the $J/\psi\pi$ invariant-mass distributions are nicely
replicated. The peak in the $D\bar{D}^{*}$ distributions is still present but
somewhat smaller. We conclude that this procedure of performing multiple joint
fits can reasonably reconcile lattice data and the main features of the
experimental data. In particular, the almost noninteracting lattice
eigenenergies can be reasonably reconciled with the experimental peaks within
the employed effective theory approach. All employed fits contain two poles
relatively close to the threshold (within
$\simeq$50\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$), and the coupling between
$D\bar{D}^{*}$ and $J/\psi\pi$ channels plays an important role.

Let us address the possibilities one could consider when trying to explain the
mentioned discrepancy between the global fit (Fit1 and Fit2) and our finite-
volume energies slightly above the $D\bar{D}^{*}$ threshold.

The EFT omits the $D^{*}\bar{D}^{*}$ channel. This, nevertheless, is not
expected to be the primary source of disagreement. It is unlikely that
including the $D^{*}\bar{D}^{*}$ channel into our Lagrangians would render
almost noninteracting energies, as observed in our simulation. In addition, we
extracted the lattice spectrum without $D^{*}\bar{D}^{*}$ and $\eta_{c}\rho$
interpolators, and it is almost unmodified. Furthermore, considering the
decays of $Z_{c}(3900)$ (near $D\bar{D}^{*}$ threshold) and $X(4020)$ (near
$D^{*}\bar{D}^{*}$ threshold) into charmed mesons, experiments do not observe
$X(4020)$ decaying into $D\bar{D}^{*}$, which is on the other hand, true for
$Z_{c}(3900)$. This suggests that the $D\bar{D}^{*}$-$D^{*}\bar{D}^{*}$
channel coupling may be insignificant.

A possible reason for the inconsistency between finite-volume energies and the
fit at higher energies could be the omission of pion-exchange contributions in
EFT, which become significant as we move to higher energies where
derivative/momentum suppression does not apply.

The origin of the discrepancy may also be inaccurate lattice energies, which
do not correspond to the continuum limit. Furthermore, the basis of
interpolating operators does not contain local diquark-antidiquark
interpolators in our simulation and in Chen  _et al._ (2019). This may be the
limiting factor, although previous works Cheung  _et al._ (2017); Padmanath
_et al._ (2015); Prelovsek  _et al._ (2015) conclude that the addition of
diquark-antidiquark interpolators does not modify the eigenenergies
significantly. Finite-volume spectra on more lattice volumes would also be
beneficial for more robust conclusions.

##  VIII Brief resume of phenomenological studies related to $1^{+-}$ and
$1^{++}$ channels

###  VIII.1 Channel with $J^{PC}=1^{+-}$ and $I=1$

In recent years, the theoretical exploration of the fascinating $Z_{c}$ states
has gained considerable attention. This includes studies on hadronic molecules
Wang  _et al._ (2013); Dong  _et al._ (2013); Wilbring  _et al._ (2013); Aceti
_et al._ (2014); Gong  _et al._ (2016); Chen  _et al._ (2023); Liu  _et al._
(2024a, b), compact tetraquark states Maiani  _et al._ (2013) and kinematic
effects from a threshold cusp Liu and Li (2013); Swanson (2015, 2016) or a
triangle singularity Liu  _et al._ (2016); Szczepaniak (2015); Chen  _et al._
(2023); Liu  _et al._ (2024a).

The JPAC collaboration performed a coupled-channel study of the $e^{+}e^{-}\to
J/\psi\pi\pi$ and $\to D\bar{D}^{*}\pi$ data from BESIII Pilloni  _et al._
(2017). They considered several amplitude parametrizations to constrain
different phenomenological models for $Z_{c}(3900)$: compact tetraquark,
hadronic molecule, virtual state and a pure kinematical enhancement. All the
models fit the data reasonably well. The best fit is obtained for a compact
tetraquark state; however, the other scenarios have not been ruled out.

An up-to-date analysis of the experimental data related to $Z_{c}(3900)$
suggests that both the molecular and nonmolecular components are comparable
Chen  _et al._ (2023). This conclusion is drawn using a generalized expression
for compositeness, which is not model-independent Matuschek  _et al._ (2021).

The internal structure of the $Z_{c}(3900)$ is also investigated employing
functional methods Hoffer  _et al._ (2024). Three potential configurations are
explored: a tight $\bar{c}c$ core surrounded by the light $\bar{q}q$ pair
($J/\psi\pi$), heavy-light meson clusters ($D\bar{D}^{*}$) and diquark-
antidiquark. Norm contributions are computed using the on-shell Bethe-Salpeter
amplitude. The results suggest that the dominant and subdominant components
are $J/\psi\pi$ and $D\bar{D}^{*}$, respectively, with a substantial
$J/\psi\pi$-$D\bar{D}^{*}$ mixing.

The vicinity of $Z_{c}$ peak to the threshold led several studies to suggest
that this peak might just be a result of a kinematical effect Swanson (2015).
However, if one wants to fit the $Z_{c}$ within a perturbative approach, the
presence of a pronounced peak calls for such a large coupling constant that
using a perturbative approach is not justified Guo  _et al._ (2015). The
resummation of the whole series leads to a near-threshold pole and qualifies
for a state Guo  _et al._ (2015).

Our findings emphasize the crucial role of the cross-channel interaction
$J/\psi\pi$-$D\bar{D}^{*}$, as illustrated in Fig. 18. The lattice studies
from the HAL QCD collaboration Ikeda  _et al._ (2016); Ikeda and for HAL QCD
Collaboration (2018) also support that this cross-channel interaction is
crucial for the existence of $Z_{c}$. A similar conclusion is drawn in Ref.
Ortega  _et al._ (2019), where the constituent quark model with coupled
channels ($J/\psi\pi$, $\eta_{c}\rho$, $D\bar{D}^{*}$, $D^{*}\bar{D}^{*}$) is
utilized. The study reveals that the $Z_{c}(3900)$ signal is due to a pole
just below the $D\bar{D}^{*}$ threshold and finds that the non-diagonal
$J/\psi\pi$-$D^{(*)}\bar{D}^{*}$ and $\eta_{c}\rho$-$D^{(*)}\bar{D}^{*}$
couplings are crucial for the $Z_{c}$ peak.

There are many works that report poles near $D\bar{D}^{*}$ threshold, either
above Aceti  _et al._ (2014); Pilloni  _et al._ (2017); Albaladejo  _et al._
(2016); Du  _et al._ (2022); Chen  _et al._ (2023) or below Guo  _et al._
(2013); Albaladejo  _et al._ (2016); He and Chen (2018); Ortega  _et al._
(2019) the threshold. Two poles near the threshold are explicitly reported in
Ref. Zhou and Xiao (2015). It utilizes a unitarized coupled-channel model and
simultaneously analyses the three invariant mass distributions of
$J/\psi\pi^{\pm}$, $(D\bar{D}^{*})^{\pm}$ and $(D^{*}\bar{D}^{*})^{\pm}$ in
the $e^{+}e^{-}$ production processes around
$4.26\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ Ablikim  _et al._ (2013, 2014a,
2014b). Their results suggest that the $Z_{c}$ signal is caused by the two
poles at $m=$3846\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$,
$\Gamma/2=$19\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ (RS II) and
$m=$3875\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$,
$\Gamma/2=$16\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ (RS III), where the one
on RS III contributes dominantly to the $Z_{c}$ peak.

###  VIII.2 Channel with $J^{PC}=1^{++}$ and $I=1$

There are not many studies on the $1^{++}$ state with isospin 1 that would be
an isospin partner of the $\chi_{c1}(3872)$. The isospin-1 state is
anticipated in the diquark-antidiquark picture Maiani  _et al._ (2005); Anwar
_et al._ (2018). The leading order of nonrelativistic EFT also predicts such a
state within the molecular scenario Ji  _et al._ (2022). A chiral EFT has
recently advanced this understanding by incorporating one-pion exchange and
three-body dynamics Zhang  _et al._ (2024). They predict the existence of a
molecular isospin-1 state related to a virtual state pole and argue that it
has not been observed in experiment since it would show up only as either a
mild cusp or dip at the threshold. Their results are similar to ours, where a
virtual state pole is found in the one-channel $D\bar{D}^{*}$ scattering (see
Sec. VI).

##  IX Elastic s-wave $J/\psi\pi$ scattering in the $1^{+-}$ channel

In this section, we examine elastic s-wave $J/\psi\pi$ scattering for
$J^{PC}=1^{+-}$. We investigate this process for two reasons. First, the
lowest threshold is $m_{J/\psi}+m_{\pi}$, and the lowest few finite-volume
states correspond to the $J/\psi\pi$ channel. Second, the strength of the
$J/\psi\pi$ interaction provides insight into whether $\pi\pi$ exchange
contributions are sufficiently attractive to render a $\bar{c}c\bar{c}c$ state
near the $J/\psi J/\psi$ threshold Dong  _et al._ (2021a), as detailed at the
end of this section. This $\pi\pi$ exchange in $J/\psi J/\psi$ scattering can
be seen in Fig. 20 where each vertex represents the $J/\psi\pi\to J/\psi\pi$
coupling.

Figure 20: A contribution to the scattering of two $J/\psi$ mesons.

Our lattice energy levels dominated by $J/\psi\pi$ interpolators reveal no
substantial interaction between $J/\psi$ and $\pi$. This is expected, as
$J/\psi$ and $\pi$ mesons do not contain any valence quarks of the same
flavor. When checking individual energy levels (see Figs. 1 and 2), one
encounters exceptions only for the few highest levels. The rest closely align
with the free energies, deviating by no more than 1$\sigma$. To make the
connection to the scattering amplitude, values for inverse $k\cot{\delta_{0}}$
are computed from our energy levels via (13) and shown in Fig. 21. Only those
below the next threshold, $\eta_{c}\rho$, are shown. They are all zero within
$1\sigma$ uncertainties, corresponding to the noninteracting scenario. This
gives us upper bounds on the scattering length $|a_{0}|$ near the $J/\psi\pi$
threshold (within the approximation $k\cot{\delta_{0}}\approx 1/a_{0}$)

$\left|a_{0}\right|\lesssim$0.1\text{\,}\mathrm{f}\mathrm{m}$$ (29)

and bound on $\left|1/(k\cot{\delta_{0}})\right|$ elsewhere
$\displaystyle\left|1/(k\cot{\delta_{0}(k)})\right|$
$\displaystyle\lesssim$0.0005\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}^{-1}$,\
E_{\textrm{cm}}-m_{J/\psi}-m_{\pi}\approx$300\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$\>,$
(30) $\displaystyle\left|1/(k\cot{\delta_{0}(k)})\right|$
$\displaystyle\lesssim$0.0015\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}^{-1}$,\
E_{\textrm{cm}}-m_{J/\psi}-m_{\pi}\approx$350\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$\>,$
which are all determined at our unphysical pion mass of
$m_{\pi}\simeq$280\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. This also
constrains the $J/\psi J/\psi$ interaction, as discussed below.

Figure 21: Inverse $k\cot{\delta_{0}}$ from s-wave $J/\psi\pi$ scattering
derived from our lattice finite-volume energies. Values around zero represent
no interaction. The energy range represents the elastic scattering region.

We also compute the $J/\psi\pi$ scattering length at physical meson masses
from the coupled-channel amplitude (20) obtained from fitting the experimental
line shapes and lattice data: the results from the different fits in Table 4
are in the range $|a_{0}|\simeq 0.01$–$0.06\text{\,}\mathrm{f}\mathrm{m}$,
which is consistent with the bound in Eq. (29) at
$m_{\pi}\simeq$280\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$.

To our knowledge, only two previous lattice determinations of the $J/\psi\pi$
scattering length are available. The quenched simulation of Yokokawa  _et al._
(2006) renders $a_{0}=$0.0119(39)\text{\,}\mathrm{f}\mathrm{m}$$, while the
dynamical simulation of Liu  _et al._ (2008); Liu (2009) gives the upper bound
$|a_{0}|\lesssim$0.01\text{\,}\mathrm{f}\mathrm{m}$$. These values are
extracted at the physical point from chiral extrapolations of three or four
different unphysical values of $m_{\pi}$. These results are significantly more
precise than ours. However, they are not directly comparable as our bound (29)
applies at $m_{\pi}\simeq$280\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. A
phenomenological evaluation from Ref. Liu  _et al._ (2013b) suggests
$|a_{0}|\lesssim$0.02\text{\,}\mathrm{f}\mathrm{m}$$. Another phenomenological
estimate that uses the central value of a quark model calculation of the
charmonium chromopolarizability Brambilla  _et al._ (2016); Dong  _et al._
(2023) and relations from Dong  _et al._ (2021a) predicts
$|a_{0}|\approx$0.04\text{\,}\mathrm{f}\mathrm{m}$$.

Let us now discuss the possible implications of $J/\psi\pi$ interaction for
the $\bar{c}c\bar{c}c$ tetraquarks. The spectroscopy of fully-charm
tetraquarks became more active in 2020 when LHCb detected the $X(6900)$ in the
$J/\psi J/\psi$ final state, supporting the $cc\bar{c}\bar{c}$ quark content
Aaij  _et al._ (2020). Subsequently, additional significant peaks were
observed, representing further states with the same proposed quark content.
There are interpretations that these states are compact tetraquarks Zhang  _et
al._ (2022); Wang  _et al._ (2022); An  _et al._ (2023); Lü  _et al._ (2020);
liu  _et al._ (2024); Yu  _et al._ (2023); Deng  _et al._ (2021); Dong and
Wang (2023); Maiani and Pilloni (2022); Ortega  _et al._ (2023); Wang  _et
al._ (2023); Anwar and Burns (2023); Nefediev (2021) or hadronic molecules
Dong  _et al._ (2021b, a); Baru  _et al._ (2022). Latter studies suggest that
the exchange of correlated light mesons, such as two pions, contributes
significantly to the attraction within a double-charmonium system. The absence
of any interaction between $J/\psi$ and $\pi$ would disfavor this scenario.

The $J/\psi\pi$ scattering length and the long-range interaction potential
between $J/\psi J/\psi$ are connected through the contact interaction vertex
$J/\psi J/\psi\pi\pi$ of the chiral effective theory (see Fig. 20). These two
different channels are related in Ref. Dong  _et al._ (2021a), where the
authors conclude that a relatively small s-wave $J/\psi\pi$ scattering length
$\mathcal{O}($0.01\text{\,}\mathrm{f}\mathrm{m}$)$ is significant enough to
explain the attraction between two $J/\psi$ mesons via the exchange of
correlated light mesons to a near-threshold state. Such small scattering
lengths are allowed, but not guaranteed, from the only available dynamical QCD
simulations that render only the upper bounds
$\left|a_{0}\right|\lesssim$0.1\text{\,}\mathrm{f}\mathrm{m}$$ (present work)
and $|a_{0}|\lesssim$0.01\text{\,}\mathrm{f}\mathrm{m}$$ Liu  _et al._ (2008);
Liu (2009). Therefore, progress on the role of two-pion exchange calls for a
more precise determination of the $J/\psi\pi$ scattering length from lattices
with dynamical quarks. This is particularly challenging due to the smallness
of this interaction.

##  X Conclusions

A charmoniumlike tetraquark $Z_{c}^{+}(3900)$ with charged decay channels and
$I(J^{PC})=1(1^{+-})$ was discovered in 2013. Its exotic quark content is
$\bar{c}c\bar{q}q$, and its peak lies just above the $D\bar{D}^{*}$ threshold.
Our investigation aims to understand whether the interaction between
$D\bar{D}^{*}$ and its coupling to $J/\psi\pi$ contributes to the existence of
$Z_{c}(3900)$ in the $1^{+-}$ sector. Furthermore, we explore whether the
$D\bar{D}^{*}$ interaction provides clues about an undiscovered near-threshold
state in the $1(1^{++})$ sector.

This work investigates the channels $\bar{c}c\bar{q}q$ with $I=1$ and
$J^{PC}=1^{+\pm}$ on the lattice. Our and other lattice finite-volume energies
for $1^{+-}$ are compared to experimental $Z_{c}(3900)$ data via an EFT. The
main challenge is that such states decay to $D\bar{D}^{*}$ as well as lower-
lying final states of a charmonium and a light meson.

In the lattice part of the study, meson-meson interpolator types (2) are
employed, and the finite-volume spectrum is extracted (see Figs. 1–5). This is
the first study considering hadronic states with $1^{+\pm}$, a non-zero total
momentum and two different lattice volumes. The determined eigenenergies are
compared to the noninteracting energies of the corresponding two-hadron
systems $\bar{c}c+\bar{q}q$ and $\bar{c}q+\bar{q}c$. The number of extracted
energy levels equals the number of expected noninteracting meson-meson levels.
All eigenenergies lie relatively close to the noninteracting meson-meson
energies, which agrees with the findings of all previous lattice studies of
these channels. Only certain eigenstates, predominantly coupled to
$D\bar{D}^{*}$, exhibit slightly negative energy shifts, which are more
pronounced in the $1^{++}$ compared to the $1^{+-}$ channel.

After extracting the finite-volume energies, our first approach is to
investigate whether single-channel $D\bar{D}^{*}$ scattering supports the
existence of $Z_{c}(3900)$ and its $1^{++}$ partner. This part of the study
employs the standard Lüscher formalism for one-channel s-wave $D\bar{D}^{*}$
scattering and the energy dependence is parametrized via the effective range
expansion. The resulting scattering amplitudes exhibit a virtual state pole
slightly below the threshold in both sectors, $1^{+-}$ and $1^{++}$. However,
we acknowledge that this approach neglects the impact of the left-hand cut and
assumes that $D\bar{D}^{*}$ does not couple to other channels, which may not
be valid assumptions. Due to these limitations, our emphasis lies in the
results discussed below.

The second procedure applies to $1^{+-}$ and involves
$J/\psi\pi$-$D\bar{D}^{*}$ coupled-channel scattering based on a contact EFT.
Considering contact interactions $D\bar{D}^{*}$-$D\bar{D}^{*}$ and
$D\bar{D}^{*}$-$J/\psi\pi$, we jointly fit the experimental $J/\psi\pi$ and
$D\bar{D}^{*}$ invariant-mass distributions from BESIII Ablikim  _et al._
(2017a, 2015b) and lattice finite-volume energies (from our work and others
Cheung  _et al._ (2017); Chen  _et al._ (2019)) with an effective field theory
approach. We perform four different fits, which differ in treating higher-
lying energy levels from our lattice study. All fits reasonably reproduce
$Z_{c}$ peak in the experimental $J/\psi\pi$ and $D\bar{D}^{*}$ line shapes.
They all reproduce also our and other lattice energy levels up to slightly
above the $D\bar{D}^{*}$ threshold. The reproduction of our higher-lying
energy levels is more favorable in fits where these levels are incorporated
with artificially reduced uncertainties. This implies that the employed EFT
model can reconcile the peaks in the experimental line shapes and the lattice
energies that lie close to noninteracting energies. All the fits yield two
poles within $\simeq$50\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ of the
$D\bar{D}^{*}$ threshold: both (RS III and RS IV) poles are possibly equally
responsible for the $Z_{c}$ peak. The pole locations of the preferred
conventional fit are shown in Fig. 16. Another common feature of all fits is
that the coupling between the $D\bar{D}^{*}$ and $J/\psi$ channels is
significant and plays a crucial role, in agreement with the findings of the
HAL QCD study Ikeda  _et al._ (2016).

Lastly, elastic s-wave $J/\psi\pi$ scattering with $1^{+-}$ is considered and
the resulting $1/(k\cot{\delta_{0}})$ is presented in Fig. 21. It is zero
within the uncertainties, which implies nearly noninteracting $J/\psi$ and
$\pi$. We, therefore, set an upper bound on the $J/\psi\pi$ scattering length
and compare it to other studies. It would be important to determine this
scattering length in the future more precisely since one could use it to
explain the binding of $J/\psi$-$J/\psi$ mesons via the exchange of correlated
light mesons to a near-threshold bound state.

The current understanding could be improved by lattice investigations of
discretization effects based on improved actions and smaller lattice spacings.
The method of optimized distillation profiles introduced in Ref. Knechtli  _et
al._ (2022) and used in Urrea-Niño  _et al._ (2023); Neuendorf  _et al._
(2024) could be applied to define operators with large overlaps onto energy
eigenstates of interest, significantly improving on standard distillation
Peardon  _et al._ (2009). Our large (statistical) uncertainties call for
ensembles with higher statistics. Since $\rho$ and $a_{0}$ resonances are
involved, one should perform three-particle scattering. To tackle this, one
has to incorporate the scattering formalism outlined, for example, in Hammer
_et al._ (2017); Mai and Döring (2017); Hansen and Sharpe (2019), where the
presence of the left-hand cut is not discussed. It is advisable to incorporate
techniques that address the left-hand cut Du  _et al._ (2023); Raposo and
Hansen (2023); Meng  _et al._ (2024) even before attempting to reduce the pion
mass. Simulating at the physical pion mass will introduce challenges related
to strongly decaying $D^{*}\to D\pi$ (consequently $Z_{c}(3900)\to
D\bar{D}\pi$). To incorporate $D\bar{D}\pi$ effects and the left-hand cut, the
method from Ref. Hansen  _et al._ (2024) can be used. Figures 12–14 with
predictions of discrete lattice spectra in a wide range of volume sizes can
provide helpful guidelines for future lattice calculations.

## Acknowledgments

We thank David R. Entem, Feng-Kun Guo, Christoph Hanhart, Luka Leskovec, Marko
Maležič, Mikhail Mikhasenko, Daniel Mohler, Raquel Molina, Alexey V. Nefediev,
Pablo G. Ortega, and Eulogio Oset for valuable discussions. We are grateful to
the members of CLS for the effort in generating the gauge field ensembles,
which form a basis for our computation. We thank Christian B. Lang and Stefano
Piemonte for contributions to the computing codes we used and the authors of
the TwoHadronsInBox package Ref. Morningstar  _et al._ (2017) for making the
TwoHadronsInBox package public. Our code implementing the distillation method
was written within the framework of the Chroma software package Edwards and
Joó (2005). The correlators were computed on two HPC systems: we gratefully
acknowledge the HPC RIVR consortium (www.hpc-rivr.si) with EuroHPC JU
(eurohpc-ju.europa.eu) for one part and the University of Regensburg for
funding this research by providing computing resources of the HPC system Vega
at the Institute of Information Science (www.izum.si) and the HPC-Cluster
Athene, respectively. M. S. acknowledges the financial support by the
Slovenian Research and Innovation Agency ARIS (Grant No. 53647). S. C.
acknowledges the support from the DFG grant SFB/TRR 55. The work of S. P. is
supported by the Slovenian Research and Innovation Agency ARIS (research core
funding No. P1-0035 and No. J1-8137) and (at the beginning of the project) DFG
grant No. SFB/TRR 55. M.P. gratefully acknowledges support from the Department
of Science and Technology, India, SERB Start-up Research Grant No.
SRG/2023/001235 and Department of Atomic Energy, India. Z. H. G. is partially
supported by the Natural Science Foundation of China (NSFC) under contracts
No. 11975090 and the Science Foundation of Hebei Normal University with
contract No. L2023B09.

##  Appendix A Masses of mesons on our lattices

The masses of single hadrons from lattice and experiment are listed in Table
6. The masses from lattice reveal a slight dependence on the lattice size
$N_{L}=24,~{}32$. Consequently, the curves representing the noninteracting
energies of both mesons with zero momentum are not precisely horizontal, as
previously elaborated upon in Sec. V.

Table 6: Masses of relevant mesons measured on both lattices and the experiment. If applicable, experimental masses of a charged and a neutral particle are shown. |  $m_{\pi^{\pm}}/m_{\pi^{0}}$ [MeV] |  $m_{D^{\pm}}/m_{D^{0}}$ [MeV] |  $m_{D^{*\pm}}/m_{D^{*0}}$ [MeV] |  $m_{J/\psi}$ [MeV] |  $m_{\eta_{c}}$ [MeV] |  $M_{\textrm{av}}$ [MeV]
---|---|---|---|---|---|---
$N_{L}=24$ |  $289(4)$ |  $1933(2)$ |  $2055(3)$ |  $3126(1)$ |  $3028(1)$ |  $3101(1)$
$N_{L}=32$ |  $280(3)$ |  $1927(1)$ |  $2049(2)$ |  $3129(1)$ |  $3030.7(3)$ |  $3104.3(4)$
exp. Workman  _et al._ (2022) |  $139.57/134.98$ |  $1869.66(5)/1864.84(5)$ |  $2010.26(5)/2006.85(5)$ |  $3096.90(1)$ |  $2983.9(4)$ |  $3068.7(1)$

##  Appendix B Lists of implemented two-hadron interpolators

All incorporated interpolators for $1^{+-}$ and $1^{++}$ are schematically
shown in Tables 7 and 8, respectively. The multiplicity of all linearly
independent interpolators is indicated by “$[\times 2]$”, “$[\times 3]$”.

Single-meson pseudoscalar and vector operators are realized using $\gamma_{5}$
and $\gamma_{i}$. In addition, specific meson-meson interpolators are
duplicated (as indicated by “$\times 2$”) by using also $\gamma_{5}\gamma_{t}$
and $\gamma_{i}\gamma_{t}$, for example

$\displaystyle O^{J/\psi(0)\pi(0)}_{S=1,\ell=0,J=1}$
$\displaystyle=\bar{c}\gamma_{z}c({\hbox{\boldmath{$0$}}})\,\bar{q}\gamma_{5}q({\hbox{\boldmath{$0$}}}),\quad\bar{c}\gamma_{z}\gamma_{t}c({\hbox{\boldmath{$0$}}})\,\bar{q}\gamma_{5}\gamma_{t}q({\hbox{\boldmath{$0$}}})\>.$

Table 7: Table of employed meson-meson interpolators (2) transforming under irreps $\Lambda^{PC}=T_{1}^{+-}$ and $\Lambda^{C}=A_{2}^{-}$ which correspond to $J^{PC}=1^{+-}$. All momenta $|{\hbox{\boldmath{$P$}}}|^{2}$, $|{\hbox{\boldmath{$p$}}}_{i}|^{2}$ are in units of $(2\pi/L)^{2}$. $|{\hbox{\boldmath{$P$}}}|^{2}=0,\ \Lambda^{PC}=T_{1}^{+-}$ | $|{\hbox{\boldmath{$P$}}}|^{2}=1,\ \Lambda^{C}=A_{2}^{-}$
---|---
|  | $J/\psi(0)\pi(0)$ | $\times 2$ | $J/\psi(1)\pi(0)$ | $\times 2$ |
|  | $J/\psi(1)\pi(1)$ | [$\times 2$] | $J/\psi(0)\pi(1)$ | $\times 2$ |
|  | $J/\psi(2)\pi(2)$ | [$\times 3$] | $J/\psi(2)\pi(1)$ | [$\times 2$] |
$N_{L}=24$ |  | $\eta_{c}(0)\rho(0)$ |  | $J/\psi(1)\pi(2)$ | [$\times 2$] |
15 interpolators | $N_{L}=32$ | $\eta_{c}(1)\rho(1)$ | [$\times 2$] | $J/\psi(4)\pi(1)$ |  |
| 21 interp. | $\bar{D}^{*}(0)D(0)$ | $\times 2$ | $\eta_{c}(1)\rho(0)$ |  | $N_{L}=24$
|  | $\bar{D}^{*}(1)D(1)$ | [$\times 2$] | $\eta_{c}(0)\rho(1)$ |  | 21 interp.
|  | $\bar{D}^{*}(0)D^{*}(0)$ |  | $\eta_{c}(2)\rho(1)$ | [$\times 2$] |
|  | $J/\psi(3)\pi(3)$ | [$\times 2$] | $\bar{D}^{*}(0)D(1)$ | $\times 2$ |
|  | $\eta_{c}(2)\rho(2)$ | [$\times 3$] | $\bar{D}^{*}(1)D(0)$ | $\times 2$ |
|  | $h_{c}(1)\pi(1)$ |  | $\bar{D}^{*}(1)D(2)$ | [$\times 2$] |
|  |  |  | $\bar{D}^{*}(2)D(1)$ | [$\times 2$] |
Table 8: Table of interpolators (2) transforming under irreps $\Lambda^{PC}=T_{1}^{++}$ and $\Lambda^{C}=A_{2}^{+}$ which correspond to $J^{PC}=1^{++}$. All momenta $|{\hbox{\boldmath{$P$}}}|^{2}$, $|{\hbox{\boldmath{$p$}}}_{i}|^{2}$ are in units of $(2\pi/L)^{2}$. $|{\hbox{\boldmath{$P$}}}|^{2}=0,\ \Lambda^{PC}=T_{1}^{++}$ | $|{\hbox{\boldmath{$P$}}}|^{2}=1,\ \Lambda^{C}=A_{2}^{+}$
---|---
$N_{L}=24$ |  | $J/\psi(0)\rho(0)$ |  | $\eta_{c}(1)a_{0}(0)$ |  |  |
$\ 5\ $ interpolators |  | $\bar{D}^{*}(0)D(0)$ | $\times 2$ | $\chi_{c0}(1)\pi(0)$ |  |  |
| $N_{L}=32$ | $\bar{D}^{*}(1)D(1)$ | [$\times 2$] | $\chi_{c0}(0)\pi(1)$ |  |  |
| 10 interp. | $J/\psi(1)\rho(1)$ | [$\times 3$] | $J/\psi(1)\rho(0)$ |  |  | $N_{L}=24$
| $\chi_{c0}(1)\pi(1)$ |  | $J/\psi(0)\rho(1)$ |  |  | 13 interp.
| $\chi_{c1}(1)\pi(1)$ |  | $\bar{D}^{*}(0)D(1)$ | $\times 2$ | $N_{L}=32$ |
|  |  | $\bar{D}^{*}(1)D(0)$ | $\times 2$ | 17 interp. |
|  | $\bar{D}^{*}(1)D(2)$ | [$\times 2$] |  |
|  | $\bar{D}^{*}(2)D(1)$ | [$\times 2$] |  |
|  | $\eta_{c}(0)a_{0}(1)$ |  |  |
|  | $\chi_{c0}(2)\pi(1)$ |  |
|  | $\chi_{c0}(4)\pi(1)$ |  |
|  | $\chi_{c1}(2)\pi(1)$ |  |

##  Appendix C Overlaps between eigenstates and operators

The correlation matrix (1) is constructed with the interpolating operators
$O_{i}$ in (2), with quantum numbers of the states we are interested in, and
can be written as

$\begin{split}C_{ij}(t)&=\langle 0|O_{i}(t)O_{j}^{\dagger}(0)|0\rangle\\\
&=\sum_{n=1}^{\infty}\langle 0|O_{i}|n\rangle\langle
n|O^{\dagger}_{j}|0\rangle e^{-E_{n}t}\\\
&\approx\sum_{n=1}^{N}Z_{i}^{(n)}Z_{j}^{(n)*}e^{-E_{n}t}\>.\end{split}$ (31)

Here, we identify an infinite tower of eigenstates of the Hamiltonian222222The
lattice eigenenergies $E_{n}^{\textrm{lat}}$ are denoted for simplicity by
$E_{n}$ in this and the following appendix., denoted as $|n\rangle$, and
normalized as $\langle m|n\rangle=\delta_{mn}$. The overlaps between the
eigenstates and the operators are expressed as:

$\begin{split}Z_{i}^{(n)}=&\langle 0|O_{i}|n\rangle\>.\end{split}$ (32)

In the final step of (31), the summation over eigenstates is truncated such
that the number of states is the same as the number of employed interpolators
($N$).

The overlaps (32) are used to identify an eigenstate $|n\rangle$, and they are
practically determined by the following equation

$Z_{i}^{(n)}(t)=e^{E_{n}t/2}\frac{\left|C_{ij}(t)u_{j}^{(n)}(t)\right|}{\left|C(t)^{\frac{1}{2}}{\hbox{\boldmath{$u$}}}^{(n)}(t)\right|}\>,$
(33)

where ${\hbox{\boldmath{$u$}}}^{(n)}$ are the eigenvectors from Eq. (8).

Since the overlaps are time-dependent, they may not be constant within the
region of interest. This challenge is particularly pronounced when dealing
with many interpolators. To address this issue, the following procedure is
adopted: The usual GEVP from Eq. (8) is solved only for $t$ up to
$t_{\mathrm{fix}}$. The eigenvectors at higher times $t>t_{\textrm{fix}}$ are
therefore not calculated. Instead, they are defined as
${\hbox{\boldmath{$u$}}}^{(n)}(t>t_{\textrm{fix}})={\hbox{\boldmath{$u$}}}^{(n)}(t_{\textrm{fix}})$.
The eigenvalues $\lambda^{(n)}(t>t_{\textrm{fix}},t_{0})$ are computed with
$C(t){\hbox{\boldmath{$u$}}}^{(n)}(t_{\textrm{fix}})=\lambda^{(n)}(t,t_{0})C(t_{0}){\hbox{\boldmath{$u$}}}^{(n)}(t_{\textrm{fix}})$
($t>t_{\textrm{fix}}$) and the overlaps $Z_{i}^{(n)}(t>t_{\textrm{fix}})$
become constant by definition. Our specific choices for $t_{\textrm{fix}}$ are
documented in Table 9. These choices are partially dictated by constraints
from mixing, described in Appendix E. Ideally, values of $t_{\textrm{fix}}$
should be consistent for all eigenstates of a given correlator. However, in
some cases, different values are necessary to disentangle some states high up
in the spectrum.

The overlaps $Z_{i}^{(n)}(t)$ are averaged over a time interval that
encompasses the plateau region of effective energies $E_{n}^{\textrm{eff}}(t)$
defined in Appendix D. The time-averaged overlaps are denoted as
$Z_{i}^{(n)}$. A valuable measure for overlaps is the normalized overlap

$\tilde{Z}_{i}^{(n)}\equiv\frac{Z_{i}^{(n)}}{\max^{N}_{m=1}Z_{i}^{(m)}}$ (34)

as it is independent of the normalization of the interpolator. An eigenstate
$|n\rangle$ is connected with the identity of an interpolator $O_{i}$ if
overlap $\tilde{Z}_{i}^{(m)}$ is the largest for eigenstate $m=n$ among all
eigenstates $m$. The choice $t_{\textrm{fix}}$ is irrelevant in many cases,
while Fig. 22 shows an example of levels where a careful choice of
$t_{\textrm{fix}}$ is valuable.

Figure 22: Examples of time-dependent overlaps $Z$ and normalized overlaps
$\tilde{Z}$ for the state $n=4$ (above) and $n=5$ (below) of the irrep
$T_{1}^{+-}$ and $N_{L}=24$. These levels highlight an example where a careful
choice $t_{\textrm{fix}}$ is valuable. In the first column, overlaps are not
fixed over time, whereas for the plot with the light yellow (green)
background, they are fixed at $t_{\textrm{fix}}=7\,(11)$. Legend represents
all the interpolators. Color and symbol schemes are similar to those from
other figures. Because of the crossover between interpolators
$D(0)\bar{D}^{*}(0)$ and $\eta_{c}(0)\rho(0)$ (first column), a different
choice of $t_{\textrm{fix}}$ manifests into different dominance. The state
$n=4$ is, e.g., dominated by $D(0)\bar{D}^{*}(0)$ when $t_{\textrm{fix}}=7$
and by $\eta_{c}(0)\rho(0)$ when $t_{\textrm{fix}}=11$, where the latter
$t_{\textrm{fix}}$ was the final choice. Note that fitted eigenenergies
$E_{n}^{\textrm{lat}}$ mostly do not depend on the choice of
$t_{\textrm{fix}}$. This is particularly true for the low-lying states.
Overlaps $Z_{i}^{(4)}(t)$ ($Z_{i}^{(5)}(t)$) are averaged over time interval
$[8,20]$ ($[10,20]$).

##  Appendix D Fitting finite-volume energy levels

The eigenvalues, $\lambda^{(n)}(t,t_{0})\propto
e^{-E_{n}^{\textrm{lat}}(t-t_{0})}$, from (8) have an asymptotically
exponential behaviour at large $t$ and give the effective energies,
$E^{(n)}_{\textrm{eff}}(t,t_{0})\equiv\frac{1}{a}\ln[\lambda^{(n)}(t,t_{0})/\lambda^{(n)}(t+a,t_{0})]$,
which equal the finite-volume eigenenergies
$E_{n}=\lim_{t\to\infty}E^{(n)}_{\textrm{eff}}(t,t_{0})$ in the plateau
region. An illustrative example of effective energies, as well as the $E_{n}$,
extracted asymptotically through single-exponential fits to
$\lambda^{(n)}(t,t_{0})$, is plotted in Fig. 23. Our choice of $t_{0}$ is
shown in Table 9. No temporal boundary effects are visible in our variational
analysis.

Figure 23: Effective energies $E^{(n)}_{\textrm{eff}}(t,t_{0}=3)$ for the irrep $T_{1}^{++}$ on the smaller lattice ($N_{L}=24$) and all eigenstates $n=1,\ldots,4$. They render eigenenergies ($E_{n}^{\textrm{lat}}$) in the plateau region shown by horizontal lines with gray bands representing their uncertainties. States dominantly couple to the interpolators shown in the legend. Dotted lines represent the noninteracting energies $E_{\textrm{ni}}^{\textrm{lat}}$. Table 9: Table with $t_{0}$ and $t_{\textrm{fix}}$ used when extracting the eigenenergies of all states below the highest states coupled to $D\bar{D}^{*}$ from the full correlation matrix. Symbol † means that those numbers are valid for all states except those listed below.  irrep |  $N_{L}$ |  states |  $t_{0}$ |  $t_{\textrm{fix}}$
---|---|---|---|---
$T_{1}^{+-}$ |  24 |  all† |  2 |  11
|  |  $J/\psi(2)\pi(2)$ s-wave |  3 |  11
|  |  $D(1)\bar{D}^{*}(1)$ s-wave |  3 |  11
|  |  $D(1)\bar{D}^{*}(1)$ d-wave |  2 |  14
|  |  $\eta_{c}(1)\rho(1)$ s-wave |  2 |  7
|  |  $\eta_{c}(1)\rho(1)$ d-wave |  2 |  16
|  |  $D^{*}(0)\bar{D}^{*}(0)$ |  3 |  7
$T_{1}^{+-}$ |  32 |  all† |  2 |  10
|  |  $D(1)\bar{D}^{*}(1)$ s-wave |  2 |  6
|  |  $D(1)\bar{D}^{*}(1)$ d-wave |  2 |  9
|  |  $\eta_{c}(1)\rho(1)$ s-wave |  2 |  9
|  |  $\eta_{c}(1)\rho(1)$ d-wave |  2 |  8
|  |  $D^{*}(0)\bar{D}^{*}(0)$ |  2 |  9
|  |  $J/\psi(3)\pi(3)$ s-wave |  2 |  9
|  |  $J/\psi(3)\pi(3)$ d-wave |  2 |  9
$T_{1}^{++}$ |  24 |  all |  3 |  16
$T_{1}^{++}$ |  32 |  all |  2 |  10
$A_{2}^{-}$ |  24 |  all |  2 |  13
$A_{2}^{+}$ |  24 |  all† |  2 |  16
|  |  $J/\psi(0)\rho(1)$ |  3 |  8
$A_{2}^{+}$ |  32 |  all† |  2 |  8
|  |  $\eta_{c}(1)a_{0}(0)$ |  3 |  10
|  |  $\chi_{c0}(2)\pi(1)$ |  2 |  10
|  |  $D(1)\bar{D}^{*}(2)$ |  3 |  10

##  Appendix E Mixing between the states

In most cases, identifying states via overlaps with specific interpolators is
clean. Exceptions do not strictly appear in cases where both states have
similar energies. For example, one would expect to encounter mixing between
$D(1)\bar{D}^{*}(1)$ and $J/\psi(3)\pi(3)$ in $T_{1}^{+-}$, $N_{L}=32$, but
there is none. A physically important example featuring relatively mixed
states is presented in the left-hand side plots in Fig. 24, which exhibits an
admixture of $D(1)\bar{D}^{*}(1)$ s-wave and $J/\psi(2)\pi(2)$ s-wave
($T_{1}^{+-}$, $N_{L}=24$). Energies from two different identifications are
shown in the right pane, and they agree within the uncertainties. Similarly,
the two interpolators $D(1)\bar{D}^{*}(0)$ and $J/\psi(1)\pi(2)$ mix in two
states of the irrep $A_{2}^{-}$ ($N_{L}=24$). These examples suggest that the
interaction between $J/\psi\pi$ and $D\bar{D}^{*}$ could be significant in
understanding experimental $Z_{c}(3900)$ peaks.

Figure 24: Two states containing significant admixtures of
$D(1)\bar{D}^{*}(1)$ s-wave, $J/\psi(2)\pi(2)$ s-wave (but also
$\eta_{c}(1)\rho(1)$ d-wave) are presented. They come from the irrep
$T_{1}^{+-}$ and $N_{L}=24$. Both panels’ color and symbol schemes are similar
to those in the other figures. Left panel: normalized overlaps for the two
states $n=7$ (above) and $n=11$ (below). Right panel: Eigenenergies from two
identifications: preferred identifications where $n\\!=\\!7$ is
$J/\psi(2)\pi(2)$ and $n\\!=\\!11$ is $D(1)\bar{D}^{*}(1)$, and alternative
identification where they are reversed. Both lead to eigenenergies that agree
within errors.

##  Appendix F More fits based on EFT

The EFT analysis of the $Z_{c}$ channel in Section VII was based on four fits
listed in Subsec. VII.2. In addition, we also performed several variations of
Fit3 with artificially reduced errors on our energy levels above
$E_{\mathrm{cm}}=$4050\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. They involve
the 4 fitted lattice points (green squares in Figs. 12, 13 and 14) with
uncertainties artificially lowered by a factor of 2, 7 and 10. Their results
show qualitatively similar behavior to Fit3. General observations are that
with more artificially enhanced precision, we get:

  * •

slightly larger $\chi^{2}/N_{\textrm{dof}}$;

  * •

better replication of our lattice energy points;

  * •

a smaller, but still present, peak of the predicted $D\bar{D}^{*}$ line
shapes;

  * •

a more significant difference in specific fitted parameters compared to Fit1,
Fit2 and Fit4.

All those fits (together with Fit3) reproduce the $J/\psi\pi$ line shapes very
similarly to other fits.

## References

  * Choi  _et al._ (2003) S. K. Choi _et al._ (Belle), Observation of a Narrow Charmoniumlike State in Exclusive ${B}^{\pm{}}\rightarrow{K}^{\pm{}}{\pi}^{+}{\pi}^{-}{J}/\psi$ Decays, Phys. Rev. Lett. 91, 262001 (2003).
  * Ablikim  _et al._ (2013) M. Ablikim _et al._ (BESIII), Observation of a Charged Charmoniumlike Structure in $e^{+}e^{-}$ → $\pi^{+}\pi^{-}$ J/$\psi$ at $\sqrt{s}$ =4.26 GeV, Phys. Rev. Lett. 110, 252001 (2013), arXiv:1303.5949 [hep-ex] .
  * Liu  _et al._ (2013a) Z. Q. Liu _et al._ (Belle), Study of $e^{+}e^{-}\to\pi^{+}\pi^{-}J/\psi$ and Observation of a Charged Charmoniumlike State at Belle, Phys. Rev. Lett. 110, 252002 (2013a), [Erratum: Phys.Rev.Lett. 111, 019901 (2013)], arXiv:1304.0121 [hep-ex] .
  * Ablikim  _et al._ (2014a) M. Ablikim _et al._ (BESIII), Observation of a charged $(D\bar{D}^{*})^{\pm}$ mass peak in $e^{+}e^{-}\to\pi D\bar{D}^{*}$ at $\sqrt{s}=$ 4.26 GeV, Phys. Rev. Lett. 112, 022001 (2014a), arXiv:1310.1163 [hep-ex] .
  * Ablikim  _et al._ (2015a) M. Ablikim _et al._ (BESIII), Observation of a Neutral Structure near the $D\bar{D}^{*}$ Mass Threshold in $e^{+}e^{-}\to(D\bar{D}^{*})^{0}\pi^{0}$ at $\sqrt{s}$ = 4.226 and 4.257 GeV, Phys. Rev. Lett. 115, 222002 (2015a), arXiv:1509.05620 [hep-ex] .
  * Ablikim  _et al._ (2017a) M. Ablikim _et al._ (BESIII Collaboration), Determination of the Spin and Parity of the ${Z}_{c}(3900)$, Phys. Rev. Lett. 119, 072001 (2017a).
  * Workman  _et al._ (2022) R. L. Workman _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * Wang  _et al._ (2013) Q. Wang, C. Hanhart, and Q. Zhao, Decoding the riddle of $Y(4260)$ and $Z_{c}(3900)$, Phys. Rev. Lett. 111, 132003 (2013), arXiv:1303.6355 [hep-ph] .
  * Dong  _et al._ (2013) Y. Dong, A. Faessler, T. Gutsche, and V. E. Lyubovitskij, Strong decays of molecular states Z${}_{c}^{+}$ and Z${}_{c}^{{}^{\prime}+}$, Phys. Rev. D 88, 014030 (2013), arXiv:1306.0824 [hep-ph] .
  * Wilbring  _et al._ (2013) E. Wilbring, H. W. Hammer, and U. G. Meißner, Electromagnetic Structure of the $Z_{c}(3900)$, Phys. Lett. B 726, 326 (2013), arXiv:1304.2882 [hep-ph] .
  * Aceti  _et al._ (2014) F. Aceti, M. Bayar, E. Oset, A. M. Torres, K. P. Khemchandani, J. M. Dias, F. S. Navarra, and M. Nielsen, Prediction of an ${I}=1$ ${D}{{\overline{D}}^{*}}$ state and relationship to the claimed ${Z}_{c}(3900)$, ${Z}_{c}(3885)$, Phys. Rev. D 90, 016003 (2014).
  * Gong  _et al._ (2016) Q.-R. Gong, Z.-H. Guo, C. Meng, G.-Y. Tang, Y.-F. Wang, and H.-Q. Zheng, $Z_{c}(3900)$ as a $D\bar{D}^{*}$ molecule from the pole counting rule, Phys. Rev. D 94, 114019 (2016), arXiv:1604.08836 [hep-ph] .
  * Chen  _et al._ (2023) Y.-H. Chen, M.-L. Du, and F.-K. Guo, Precise determination of the pole position of the exotic $Z_{c}(3900)$, (2023), arXiv:2310.15965 [hep-ph] .
  * Liu  _et al._ (2024a) M.-Z. Liu, X.-Z. Ling, and L.-S. Geng, Productions of $X(3872)$/$Z_{c}(3900)$ and $X_{2}(4013)$/$Z_{c}(4020)$ in $Y(4220)$ and $Y(4360)$ decays, (2024a), arXiv:2404.07681 [hep-ph] .
  * Liu  _et al._ (2024b) Z.-W. Liu, J.-X. Lu, M.-Z. Liu, and L.-S. Geng, Femtoscopy can tell whether $Z_{c}(3900)$ and $Z_{cs}(3985)$ are resonances or virtual states, (2024b), arXiv:2404.18607 [hep-ph] .
  * Maiani  _et al._ (2013) L. Maiani, V. Riquer, R. Faccini, F. Piccinini, A. Pilloni, and A. D. Polosa, A $J^{PG}=1^{++}$ Charged Resonance in the $Y(4260)\to\pi^{+}\pi^{-}J/\psi$ Decay?, Phys. Rev. D 87, 111102 (2013), arXiv:1303.6857 [hep-ph] .
  * Liu and Li (2013) X.-H. Liu and G. Li, Exploring the threshold behavior and implications on the nature of Y(4260) and Zc(3900), Phys. Rev. D 88, 014013 (2013), arXiv:1306.1384 [hep-ph] .
  * Swanson (2015) E. S. Swanson, $Z_{b}$ and $Z_{c}$ Exotic States as Coupled Channel Cusps, Phys. Rev. D 91, 034009 (2015), arXiv:1409.3291 [hep-ph] .
  * Swanson (2016) E. S. Swanson, Cusps and Exotic Charmonia, Int. J. Mod. Phys. E 25, 1642010 (2016), arXiv:1504.07952 [hep-ph] .
  * Liu  _et al._ (2016) X.-H. Liu, M. Oka, and Q. Zhao, Searching for observable effects induced by anomalous triangle singularities, Phys. Lett. B 753, 297 (2016), arXiv:1507.01674 [hep-ph] .
  * Szczepaniak (2015) A. P. Szczepaniak, Triangle Singularities and XYZ Quarkonium Peaks, Phys. Lett. B 747, 410 (2015), arXiv:1501.01691 [hep-ph] .
  * Pilloni  _et al._ (2017) A. Pilloni, C. Fernandez-Ramirez, A. Jackura, V. Mathieu, M. Mikhasenko, J. Nys, and A. P. Szczepaniak (JPAC), Amplitude analysis and the nature of the Zc(3900), Phys. Lett. B 772, 200 (2017), arXiv:1612.06490 [hep-ph] .
  * Hoffer  _et al._ (2024) J. Hoffer, G. Eichmann, and C. S. Fischer, Hidden-flavor four-quark states in the charm and bottom region, Phys. Rev. D 109, 074025 (2024), arXiv:2402.12830 [hep-ph] .
  * Ikeda  _et al._ (2016) Y. Ikeda, S. Aoki, T. Doi, S. Gongyo, T. Hatsuda, T. Inoue, T. Iritani, N. Ishii, K. Murano, and K. Sasaki (HAL QCD), Fate of the Tetraquark Candidate ${Z}_{c}(3900)$ from Lattice QCD, Phys. Rev. Lett. 117, 242001 (2016).
  * Ikeda and for HAL QCD Collaboration (2018) Y. Ikeda and for HAL QCD Collaboration (HAL QCD), The tetraquark candidate from dynamical lattice QCD simulations, Journal of Physics G: Nuclear and Particle Physics 45, 024002 (2018).
  * Ortega  _et al._ (2019) P. G. Ortega, J. Segovia, D. R. Entem, and F. Fernández, The ${Z}_{c}$ structures in a coupled-channels model, The European Physical Journal C 79, 78 (2019).
  * Prelovsek and Leskovec (2013a) S. Prelovsek and L. Leskovec, Search for $Z^{+}_{c}$(3900) in the $1^{+-}$ Channel on the Lattice, Phys. Lett. B 727, 172 (2013a), arXiv:1308.2097 [hep-lat] .
  * Prelovsek  _et al._ (2015) S. Prelovsek, C. B. Lang, L. Leskovec, and D. Mohler, Study of the ${Z}_{c}^{+}$ channel using lattice QCD, Phys. Rev. D 91, 014504 (2015).
  * Chen  _et al._ (2014) Y. Chen _et al._ (CLQCD), Low-energy scattering of the $({D}{{\overline{D}}^{*}}{)}^{\pm{}}$ system and the resonance-like structure ${Z}_{c}(3900)$, Phys. Rev. D 89, 094506 (2014).