Explaining the Many Threshold Structures in the Heavy-Quark Hadron Spectrum

Introduction.—Quantum chromodynamics, the fundamental theory of the strong interaction, is nonperturbative at low energies. As a result, the low-energy strong interaction is notoriously difficult to be tackled with, and so far the mechanism for the color confinement, i.e., all quarks and gluons are confined inside color neutral hadrons, remains the most challenging problem in the standard model. As a manifestation of color confinement, how the spectrum of hadrons emerges from the underlying strong dynamics remains mysterious.

One of the most puzzling issues in hadron spectroscopy is how the plethora of resonant structures observed since 2003 can be understood. Most of them were observed in the invariant mass of hadrons containing heavy quarks. The meson(-like) structures in the heavy quarkonium mass region are called the $XYZ$ states to emphasize that the internal structure and why they appear in the mass spectrum have not been understood (see Refs. Chen et al. (2016); Esposito et al. (2017); Hosaka et al. (2016); Richard (2016); Lebed et al. (2017); Guo et al. (2018); Ali et al. (2017); Olsen et al. (2018); Altmannshofer et al. (2019); Cerri et al. (2019); Liu et al. (2019); Brambilla et al. (2020); Guo et al. (2020); Yang et al. (2020) for recent reviews).

One salient feature of the new resonant hadron structures is that many of them have masses around the thresholds of a pair of hadrons. To name a few famous examples, let us mention the $X(3872)$ Choi et al. (2003) and $Z_{c}(3900)^{\pm}$ Ablikim et al. (2013a); Liu et al. (2013); Ablikim et al. (2014a) around the $D\bar{D}^{*}$ threshold, the $Z_{c}(4020)^{\pm}$ Ablikim et al. (2013b, 2014b) near the $D^{*}\bar{D}^{*}$ threshold, the $Z_{b}(10610)^{\pm}$ and $Z_{b}(10650)^{\pm}$ Bondar et al. (2012); Garmash et al. (2016) near the $B\bar{B}^{*}$ and $B^{*}\bar{B}^{*}$ thresholds, the $Z_{cs}(3985)^{-}$ Ablikim et al. (2021) near the $\bar{D}_{s}D^{*}$ and $\bar{D}_{s}^{*}D$ thresholds, the $P_{c}$ states Aaij et al. (2019) near the $\bar{D}^{(*)}\Sigma_{c}$ thresholds, and the dip and peak structures in the double $J/\psi$ spectrum Aaij et al. (2020) near the $J/\psi\psi(2S)$ and $J/\psi\psi(3770)$ thresholds, respectively. These structures were suggested to be associated with threshold cusps Bugg (2004, 2011); Chen and Liu (2011); Chen et al. (2013); Swanson (2015, 2016); Ikeda et al. (2016); Kuang et al. (2020); Sombillo et al. (2020); Dong et al. (2021); Wang et al. (2020); Karliner and Rosner (2020); Cao et al. (2020) (for a review, see Ref. Guo et al. (2020)). Notice that a pronounced threshold cusp requires the existence of a nearby pole Bugg (2008); Guo et al. (2015). However, so far it is not known when and at what thresholds, a nontrivial structure will show up. In this Letter, we aim to answer this question. We will discuss the condition to have a peak at threshold, and show that it is natural to expect structures in the final state of a heavy quarkonium and light hadrons at the threshold of a pair of open-heavy-flavor hadrons that have attractive interaction at threshold.

Single channel.—Let us start with the effective range expansion for the $S$-wave amplitude of a two-body scattering (see Refs. Guo et al. (2018); Brambilla et al. (2020) for related discussions)

where $a_{0}$ and $r_{0}$ are the $S$-wave scattering length and effective range, $k$ is the magnitude of the c.m. momentum, and $\beta$ is some hard scale of the order of the inverse of force range. Here the sign convention for the scattering length is such that it is negative for repulsive and positive for attractive interaction in the absence of a bound state. In the near-threshold region, we have the nonrelativistic expression for the momentum $k=\sqrt{2\mu E}$ with $\mu$ the reduced mass and $E$ the energy relative to the two-body threshold. One sees that the amplitude possesses a square-root branch point at $E=0$. As a result, there is a cusp exactly at the threshold in the modulus of the amplitude as a function of energy.

If we focus on the region in the immediate vicinity of the threshold, we can neglect the effective range term. Writing the amplitude as a function of $E$,

we get

It is easy to see that for positive $a_{0}>0$ (attractive interaction but not strong enough to form a bound state) this distribution is maximal, and thus has a cuspy peak, at the threshold $E=0$. Notice that in this case, there is a virtual state pole in the second Riemann sheet of the complex energy plane at $E_{\text{virtual}}=-1/(2\mu a_{0}^{2})$.

In order to check what determines the shape of the cusp, we may compute the half-maximum width of $|f_{0}(E)|$. It is ready to be obtained from the solutions of $|f_{0}(E)|=|f_{0}(0)|/2$, $E_{+}=\frac{3}{2\mu a_{0}^{2}}$ and $E_{-}=-\frac{1}{2\mu a_{0}^{2}}$, and the half-maximum width is

One sees that the cusp in the energy distribution is narrower for a larger scattering length (thus stronger attraction) and also a larger reduced mass. The former feature has been used to precisely measure the $S$-wave $\pi\pi$ scattering lengths Pislak et al. (2001, 2003); Batley et al. (2008). The latter would imply that for a fixed scattering length, the signal of the threshold cusp becomes more and more evident when we increase the reduced mass, and it is natural to observe them in the heavy-hidden-flavor sector.

For a negative $a_{0}$, the distribution above threshold decreases monotonically in exactly the same way as the case of a positive $a_{0}$; below threshold, the maximum is located at the pole in the first Riemann sheet, $E_{\text{bound}}=-1/(2\mu a_{0}^{2})$. For strong attraction, the pole is close to the threshold, and leads to a near-threshold peak. For repulsive interaction, the value of $|a_{0}|$ is small compared to the range of forces; the pole is far away beyond the applicability region of the scattering length approximation (and thus is unphysical), and no nontrivial near-threshold structure exists.

Although we showed that the energy distribution of $|f_{0}(E)|$ has a maximum at the threshold for a positive $a_{0}$, the part below threshold needs to be observed in the final state of a lower channel. Thus, we need to consider a coupled-channel problem (see Ref. Baru et al. (2005) for an analysis of the Flatté parametrization) to fully exhibit the threshold cusp structure.

Coupled channels.—Let us consider the energy region around the highest threshold in a coupled-channel system. Suppose that all open channels with lower thresholds are relatively far away so that the momentum variation for these channels around the highest threshold is smooth. In the spirit of the optical potential we may parametrize the scattering amplitude for the channel of interest in terms of a complex scattering length, i.e., Eq. (2) with $a_{0}$ taking a complex value. Unitarity requires $a_{0}$ to satisfy (see below)

With a complex $a_{0}$, $|f_{0}(E)|^{2}$ becomes

One sees that the line shape decreases monotonically above threshold because of the unitarity constraint, Eq. (5), and it also decreases below threshold if the real part of $1/a_{0}$ is positive, or $\text{Re}\,(a_{0})>0$. The half-maximum width of $|f_{0}(E)|$ in this case is

where the sum runs over $x=\text{Im}\,(1/a_{0})$ and $\text{Re}\,(1/a_{0})$. Again the cusp width is inversely proportional to the reduced mass.

To be more specific, let us consider a two-channel (denoted by channel 1 and channel 2) problem and focus on the energy region around the threshold of channel 2, which is the higher one, such that $|E|<\Delta$ with $E$ the c.m. energy relative to the higher threshold and $\Delta$ the difference between the two thresholds. One is ready to construct a nonrelativistic effective field theory (NREFT). The c.m. momentum of channel 2 is given by

The c.m.momentum of channel 1 can be expanded in power series of $E$, corresponding to even powers of $k_{2}$, as

where $m_{i,1}$ and $m_{i,2}$ are the masses of particles in channel-$i$, $\Sigma_{1}$ and $\Sigma_{2}$ are the thresholds of the lower and higher channels, respectively, and $\sqrt{s}=\Sigma_{2}+E$. At leading order (LO) of the $E$ expansion, the Bethe-Salpeter integral equation becomes an algebraic equation, and the $T$ matrix can be written as

where $G^{\Lambda}(E)$ is a diagonal matrix containing the Green’s functions for both channels. Around the channel-2 threshold, we treat the propagators of particles in this channel nonrelativistically, while those for the lower channel may be kept relativistically. That is,

with $P^{2}=s$ and $\mu_{2}$ the reduced mass of channel 2. Both loop integrals are ultraviolet divergent, and need to be regularized. Here, the nonrelativistic loop has been regularized using a hard cutoff $\Lambda_{2}>k_{2}$, for which higher order terms of $\mathcal{O}(k^{2}/\Lambda_{2}^{2})$ will be neglected, while a different regularization method may be applied to the relativistic one. For the latter, we do not show the explicit expression, but notice that the real part, $R(\Lambda_{1})$, can be expanded in Taylor series of $E$. Although $k_{1}$ can also be expanded as in Eq. (11), we keep it explicitly as it is regulator independent.

At LO, $V^{\Lambda}$ is a symmetric constant matrix, so does its inverse. The regulator-dependent terms in $G_{1}^{\Lambda}$ and $G_{2}^{\Lambda}$ can be renormalized by the diagonal matrix elements of $1/V^{\Lambda}$. Hence, we can write the $T$ matrix, which is cutoff independent, as

with $\text{det}=\left(\frac{1}{a_{11}}-i\,k_{1}\right)\left(\frac{1}{a_{22}}+\sqrt{-2\mu_{2}E-i\epsilon}\right)-\frac{1}{a_{12}^{2}}$. It is rather similar to the coupled-channel NREFT constructed in Ref. Cohen et al. (2004) (see also Ref. Braaten and Kusunoki (2005)), where both channels are treated nonrelativistically. When the channel-coupling parameter $1/a_{12}$ vanishes, the system is reduced to two single channels.

Denoting the scattering length for the higher channel including the effects from the open channel as $a_{22,\text{eff}}$, one has

The statement in Eq. (5) is proven.

For the production of the two particles in channel 1 in the energy region close to the threshold of channel 2 (see Fig. 1), the amplitude may be written as

where $V_{11,12}^{\Lambda}$ are the elements of the $V^{\Lambda}$ matrix in Eq. (12), and the short-distance factors $P_{1,2}$ are $\Lambda$-independent, analytic in $E$ and can be taken as constants at LO in the $E$ expansion. For the second term with $G_{2}^{\Lambda}$, we have used $\Lambda_{2}>k_{2}$, and kept only the cutoff term in $G_{2}^{\Lambda}$; it gets multiplicatively renormalized by the prefactor $P_{1}^{\Lambda}(V_{11}^{\Lambda})^{-1}V_{12}^{\Lambda}+P_{2}^{\Lambda}$ which should scale as $1/\Lambda_{2}$ Braaten and Kusunoki (2005). The final expression is cutoff independent.

From analyzing the $E$ dependence of the modulus squared of the amplitude in Eq. (17), we can deduce the behavior of the invariant mass distribution of particles in the lower channel around the higher threshold. We discuss two cases in the following. For illustration, we will show some line shapes for a system with the lower channel being the $J/\psi\pi^{-}$ and the higher one being the $D^{0}D^{*-}$. [Note that we only use their masses. For analyses of the system related to the $Z_{c}(3900)$ considering more complicated dynamics, see Refs. Wang et al. (2013); Albaladejo et al. (2016); Pilloni et al. (2017); Gong et al. (2018); Danilkin et al. (2020).]

Case 1: The production rate of channel 2 is much larger than that of channel 1, so that the production of channel 1 mainly proceeds through the last diagram in Fig. 1, i.e., channel 2 is the driving channel. The event distribution is proportional to $|T_{21}(E)|$. We have

Since the only nontrivial energy dependence comes from the nonanalytic, second piece in the square brackets in Eq. (18), the analysis of Eq. (8) with a complex scattering length applies. Thus, the distribution maximizes at the higher threshold in its vicinity if the real part of $a_{22,\text{eff}}$ is positive. This happens if the interaction for channel 2—which contains two parts: direct interaction in channel 2 and corrections from the channel coupling—is attractive, but not sufficient to produce a bound state. The red solid and blue dashed curves in Fig. 2 (a) illustrate such a situation. The corresponding poles in the complex $k_{2}$ plane are at $(0.04-i\,0.08)$ GeV and $(0.37-i\,0.08)$ GeV, respectively. When the pole is closer to the threshold of channel 2, the cusp is sharper.

It could happen that $a_{22}>0$ is large and $a_{11}$ is positive so that the channel coupling may render $\text{Re}\,(a_{22,\text{eff}})$ negative, or $a_{22}<0$ and $\text{Re}\,(a_{22,\text{eff}})$ remains negative after the channel coupling. In that case, the threshold would not be a local maximum, and there would be a peak below the higher threshold due to a below-threshold pole. An illustration is shown as the green dotted curve in Fig. 2 (a). The pole corresponding to these parameters is at $(-0.09-i\,0.08)$ GeV in the complex $k_{2}$ plane.

The effective scattering length $a_{22,\text{eff}}$ can be extracted using Eq. (18) from a precise measurement of the threshold structure in this case.

Case 2: The production rate of channel 1 is much larger than that of channel 2, so that the production of channel 1 mainly proceeds through the first two diagrams in Fig. 1. Then, the energy dependence of the production amplitude is dominated by that of $T_{11}(E)$. We have

One sees that if the channel coupling is weak (i.e., $a_{12}^{2}$ is large), the energy dependence around the higher threshold will be weak. However, if the channel coupling is strong, there will be nontrivial energy dependence. For $E<0$,

Thus, for $a_{22}>0$ and $1/a_{11}\leq a_{22}/a_{12}^{2}$, $|T_{11}(E)|$ increases monotonically when $E$ decreases from 0 (the red solid curve in Fig. 2 (b) presents an example). For $E\geq 0$,

The second term in the square brackets maximizes at $E=1/(2\mu_{2}a_{22}^{2})$, while the first decreases for increasing $E$ when $a_{22}>0$ and $1/a_{11}\leq a_{22}/a_{12}^{2}$. Thus, whether $|T_{11}(E)|$ increases or decreases above the threshold of channel 2 depends on the competition between the two terms, and there must be a dip for a large $|a_{22}|$ because $T_{11}$ has a zero at $\sqrt{2\mu_{2}E}=-i/a_{22}$, see Eq. (19).

In Fig. 2 (b), we show a few typical curves for the threshold behavior of $|T_{11}|^{2}$. They are clearly different from those of $|T_{21}|^{2}$ shown in the left panel of the same figure. Although there is not necessarily a visible cusp or peak, there must be a dramatic change in the near-threshold region if there is a nearby pole in the amplitude, as can be seen from the red solid and green dotted curves. It is interesting to notice that the near-threshold behavior can be a sudden but continuous decrease around the threshold in the line shape. Such a behavior is present around the $p\bar{p}$ threshold in the $\eta^{\prime}\pi^{+}\pi^{-}$ invariant mass distribution measured by BESIII Ablikim et al. (2016).

There is an important implication to the production mechanism for a near-threshold state which is rooted in channel 2 (large $a_{22}$): if it is observed as a peak in the final state of a lower channel, the driving production channel should be channel 2. Let us consider the $f_{0}(980)$ and the $K\bar{K}$ threshold as an example. In the $J/\psi\to\phi\pi^{+}\pi^{-}$ and $J/\psi\to\omega\pi^{+}\pi^{-}$ processes, the $f_{0}(980)$ should be mainly from the $s\bar{s}$ and $(u\bar{u}+d\bar{d})$ sources, respectively, and thus the $K\bar{K}$ (channel-2) and $\pi\pi$ (channel-1) as the corresponding meson channels. The driving component for the production should be $T_{21}$ and $T_{11}$, respectively. Consequently, around the $K\bar{K}$ threshold, there is a narrow peak in the $\pi^{+}\pi^{-}$ distribution for the $J/\psi\to\phi\pi^{+}\pi^{-}$ Ablikim et al. (2005) and a dip around the $K\bar{K}$ threshold for the $J/\psi\to\omega\pi^{+}\pi^{-}$ Ablikim et al. (2004). The analysis further supports that the driving channels for producing the $P_{c}$ peaks in $\Lambda_{b}$ decays are $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ Liu et al. (2016); Du et al. (2020), instead of $J/\psi p$ Roca et al. (2015) or $\Lambda_{c}\bar{D}^{(*)}$ Burns and Swanson (2019).

More complex line shapes can be produced by the interference between the two terms in Eq. (17), i.e., between case 1 and case 2, if they have comparable strengths.

Many of the near-threshold structures have been measured in final states of the threshold channels, such as the $Z_{c}(3900)$ in $D\bar{D}^{*}$ Ablikim et al. (2014a), $Z_{c}(4020)$ in $D^{*}\bar{D}^{*}$ Ablikim et al. (2014b), $Z_{cs}(3985)$ in $\bar{D}_{s}D^{*}+\bar{D}_{s}^{*}D$ Ablikim et al. (2021), $Z_{b}(10610)$ in $B\bar{B}^{*}$ Garmash et al. (2016). The production amplitude should be

with $P_{1,2}$ the short-distance parts. The event distribution in channel 2 is proportional to $k_{2}|P_{2,\text{tot}}(E)|^{2}$. Since $1/a_{11}-i\,k_{1}$ is smooth around the threshold of channel 2, the $E$ dependence of $|T_{22}(E)|$ should be almost the same as that of $T_{12}(E)=T_{21}(E)$, see Eq. (15), and thus the analysis of case 1 above applies. In the final state of channel 2, only the $E\geq 0$ piece can show up, which means that there should always be a near-threshold enhancement no matter the sign of $\text{Re}\,(a_{22,\text{eff}})$, see Eq. (8). However, the $k_{2}$ factor from the phase space would suppress the threshold behavior if the pole is not sufficiently close to the threshold. This, as well as the analysis below Eq. (18), supports the conclusion in Ref. Guo et al. (2015) in an analysis of the $Z_{c}(3900)$ that a near-threshold prominent cusp implies a near-threshold pole.

Discussions and conclusion.—From the discussions above, one expects that there should always be a peak around the threshold for a two-body system with attractive interaction in the invariant mass distribution of a lower channel (channel-1), if the production proceeds mainly through the channel with the relevant threshold, i.e., channel 2 above. This should apply to all hadron pairs with one containing a heavy quark and the other containing a heavy antiquark. Because for the heavy quark-antiquark pair to form a heavy quarkonium, the relative momentum should be small, and it should only be a small part of the whole phase space, the production of a pair of open-heavy-flavor hadrons should be much easier than that for a heavy quarkonium plus light hadrons. Consequently, it is natural to expect that the mechanism shown as the last diagram in Fig. 1 should be significant for the production of a heavy quarkonium plus light hadrons when open-flavor thresholds are open, and there should be a near-threshold peak if the interaction in the open-flavor channel is attractive. The peak is a threshold cusp if $\text{Re}\,(a_{22,\text{eff}})$ is positive, and is a peak just below threshold if the attraction is strong enough to produce a below-threshold bound state such that $\text{Re}\,(a_{22,\text{eff}})$ becomes negative and its absolute value is large. It is interesting to notice that with the same effective scattering length, the larger the reduced mass is, the sharper the peak would be.

Notice that the threshold behavior gets smeared if at least one of the relevant intermediate particles carries a finite width. Hence, it should be prominent only for those with negligible or tiny widths. The existence of a triangle singularity in special cases, when the kinematics is such that the singularity is close to the physical region, can lead to additional nontrivial structure at or above threshold (for a detailed discussion, see the review Guo et al. (2020)).

In Table 1, we list heavy-antiheavy hadron pairs (taking the charmed ones for example) that are expected to have attractive interaction at threshold. The conclusion is based on the one-boson exchange model with Lagrangians constructed considering heavy quark spin symmetry Yan et al. (1992); Casalbuoni et al. (1997); Liu and Oka (2012). The pseudoscalar exchanges in the chiral limit always yield potentials proportional to the square of transferred momentum, ${\bm{q}}^{2}$, which vanish at threshold. Hence, only vector-meson exchanges are considered, which is analogous to the vector-meson dominance contribution to the low-energy constants in chiral perturbation theory Ecker et al. (1989) (see the Supplemental Material for a list of the potentials). Calculations of heavy-antiheavy hadron interactions using lattice QCD are indispensable to reach model-independent results on which hadron pairs are attractive, and thus to understanding whether there should be a nontrivial structure around a given threshold.

It is expected that near-threshold peaks for these hadron pairs will show up in final states of a heavy quarkonium plus light hadrons. One also notices that the Born term contribution from the light-boson exchange to the scattering length scales as the heavy quark mass $m_{Q}$, and one would expect it stronger in the bottom-antibottom than in the charm-anticharm sector.

The predictions made here are ready to be tested, and the threshold structures will play an important role in revealing the interactions between a pair of heavy hadrons, and thus to the mysterious hadron structures with a pair of heavy quark and antiquark. At last, the NREFT analysis of the threshold behavior is general enough and may find its applications in fields other than hadron physics.