Over the past two decades, there has been a significant increase in the discovery of non-conventional hadrons through high energy experiments worldwide. A majority of these hadrons have been observed in the invariant mass distributions of hadrons that possess heavy quarks. The presence of meson-like structures within the heavy quarkonium mass region has led to the designation of these states as the XYZ states, indicating that their internal structure and the reasons for their appearance in the mass spectrum remain elusive and not fully understood. In recent years, other varieties of exotic states have been observed in experiments, including hidden-charm pentaquark states, double-charm tetraquark states, fully charmed tetraquark states and so on. These observations pose additional challenges, but they also create more opportunities for a better understanding of the hadron spectrum and low-energy Quantum Chromodynamics. (see Refs. [1−17] for recent reviews).

One prominent characteristic of the newly discovered resonant hadron structures is that many of them have masses near the thresholds of two hadrons. However, it remains unclear when and at which thresholds nontrivial structures will emerge. In our previous works^[18−20], our objective was to address this question. We examined the conditions necessary for a peak to occur at the threshold and demonstrated that it is reasonable to anticipate the presence of structures in the final state of a heavy quarkonium and light hadrons when they approach the threshold of a pair of open-heavy-flavor hadrons that exhibit an attractive interaction. Subsequently, we employed a one-boson-exchange model, driven by light vector meson exchange, to estimate the interaction between a pair of charmed-anticharmed hadrons and a pair of charmed-charmed hadrons, thus facilitating the search for potential molecular states within these systems.

It is a well-known fact that the unitarity of the S matrix necessitates the two-body threshold to be a square-root branch point of the scattering amplitude. Consequently, the square modulus of an S-wave amplitude will possess a cusp at the threshold, and the invariant mass distribution may exhibit a nontrivial structure. The detailed structure of the cusp, a peak, a dip, or buried in the background, is determined by the interactions between the two particles near the threshold.

To be more precise, let us consider the S-wave scattering of a two-channel (denoted by channel-1 and channel-2 with the latter having a higher threshold) case, where the scattering amplitude near the higher threshold can be cast into^[18]

with

$ a_{ij} $ the scattering length, $ \mu_2 $ the reduced mass of channel-2, $ \varSigma_i $ the threshold of channel-$ i $ and $ E=\sqrt{s}-\varSigma_2 $ the c.m. energy relative to the second threshold. Near the second threshold, the c.m. momentum of channel-2 is given by $ k_2 = \sqrt{2\mu_2 E} $ and the c.m.momentum of channel-1 can be approximated as a constant at the leading order of E. Note that we have ignored the term of $ \mathcal O(k_2^2) $, which may lead to standard Breit-Wigner-like resonance but should be negligible for the lineshape near threshold where E is small.

In the production of two particles in channel-1, with energies close to the threshold of channel-2, the amplitude can be described as

where $ P_1 $ and $ P_2 $ account for the renormalized short-distance production vertex and are constants at leading order. By analyzing the energy dependence of the modulus squared of the amplitude in Eq. (3), we can understand the behavior of the invariant mass distribution of particles in the lower channel around the higher threshold. Two cases are discussed in the following, corresponding to the two subplots in Fig. 1.

Figure 1.  Illustration of threshold behaviors. Here we use the masses of the $ \pi^- $ and $ J/\psi $ for channel-1 and those of the $ D^0 $ and $ D^{*-} $ for channel-2, and the values of used $ a_{ij} $ parameters are given in the legend. (a) Line shapes of $ |T_{21}|^2 $, which are normalized at the threshold of channel-2, i.e., $ E=0 $; (b) line shapes of $ |T_{11}|^2 $, which are normalized at $ E=-0.02 $. There is a near-threshold bound or virtual state when $ a_{22}=-3 \; {\mathrm{GeV}}$ fm or 3 fm and the attraction is much weaker when a₂₂ = 0.5 fm.

In case-1, where the production rate of channel-2 is significantly higher than that of channel-1, the event distribution is approximately proportional to $ |T_{21}(E)|^2 $. In this scenario, the near-threshold behavior of the lineshape is comparable to that of the single-channel case. If the interaction is attractive but not strong enough to form a bound state, a maximum of the amplitude square will always appear at the threshold. However, if a bound state is formed, the peak will be located at the pole position below the threshold, as shown in Fig. 1(a).

In case-2, where the production rate of channel-1 is much larger than that of channel-2, the energy dependence of the production amplitude is dominated by that of $ T_{11}(E) $. In this circumstance, instead of a peak, a dip near the threshold may emerge, as shown in Fig. 1(b). There is an essential implication for the production mechanism of a near-threshold state that stems from channel-2: if it is observed as a peak in the final state of a lower channel, then the driving production channel should be channel-2.

To illustrate this concept, we consider the $ f_0(980) $ near the $ K\bar{K} $ threshold. The $ f_0(980) $ should be primarily from the $ s\bar{s} $ and $ (u\bar{u}+d\bar{d}) $ sources in the $ J/\psi\to \phi\pi^+\pi^- $ and $ J/\psi\to \omega\pi^+\pi^- $ processes, respectively, and the corresponding meson channels are $ K\bar{K} $ (channel-2) and $ \pi\pi $ (channel-1). The driving component for the production should be $ T_{21} $ and $ T_{11} $, respectively. Consequently, around the $ K\bar{K} $ threshold, there are a narrow peak in the $ \pi^+\pi^- $ distribution for the $ J/\psi\to \phi\pi^+\pi^- $^[21] and a dip around the $ K\bar{K} $ threshold for the $ J/\psi\to \omega\pi^+\pi^- $^[22]. It is similar in the $ P_c $ case. The peaks in the $ J/\psi p $ invariant mass distribution near the $ \varSigma_c^{(*)}\bar D^{(*)} $ thresholds indicates that the driving channels for producing the $ P_c $ peaks are $ \varSigma_c^{(*)}\bar D^{(*)} $^[23–24], instead of $ J/\psi p $^[25] or $ \varLambda_c\bar D^{(*)} $^[26].

If the two terms in Eq. (3) have comparable strengths, the interference between them can produce more complex line shapes.

Based on the previous discussions, it is expected that the invariant mass distribution of a lower channel (channel-1) for a two-body system with attractive interaction should always have a peak around the threshold, if the production mainly proceeds through the channel with the relevant threshold (channel-2) above. It should apply to all hadron pairs with one containing a heavy quark and the other containing a heavy antiquark. Due to the smaller phase space of heavy quarkonium plus light hadrons, the production of a pair of open-heavy-flavor hadrons should be relatively easy, which aligns with case 1 discussed earlier. Consequently, if the interaction in the open-flavor channel is attractive, it is natural to expect a near-threshold peak. The peak would be a threshold cusp in the absence of a bound state and a peak just below threshold if the attraction is strong enough to produce a below-threshold bound state.

Table 1 lists heavy-antiheavy hadron pairs (using the example of charmed hadrons) that are expected to have attractive at threshold via vector meson exchange. This conclusion is based on the one-boson exchange model with Lagrangians constructed by considering the heavy quark spin symmetry^[27−29]. It is anticipated that near-threshold peaks for these hadron pairs will appear in the final states of heavy quarkonium plus light hadrons. Also, note that the Born term contribution from light-boson exchange to the scattering length scales as the heavy quark mass $ m_Q $ and is expected to be stronger in the bottom-antibottom sector than in the charm-anticharm sector.

The interactions at low energy between two hadrons can be described by effective field theories, with the associated low energy constants (LECs) therein in principle should be determined by experimental data. When experimental information is lacking, these LECs can be estimated through phenomenological approaches. At the leading order, the interaction potential between heavy and anti-heavy hadrons remains constant. Our study delved into the potential of heavy-antiheavy systems by exploring the resonance saturation of the contact interaction through the exchange of one-boson (specifically light pseudoscalar and vector mesons), in conjunction with heavy quark spin symmetry, chiral symmetry, and SU(3) flavor symmetry. The resonance saturation technique is recognized for effectively approximating the LECs in higher-order Lagrangians of chiral perturbation theory^[31–32]. It was observed that the values of the LECs are predominantly influenced by the exchange of vector mesons around the mass of the $\rho$-meson when vector meson exchange is permissible. Owing to the symmetry of heavy quark flavors, the potentials between pairs of bottomed hadrons mirror those between charmed hadrons. Therefore, we opted to illustrate our points using charmed hadrons as an example.

The constant potentials acquired are subsequently employed in solving the Bethe-Salpeter equation given by

where G denotes the one-loop two-body propagator that undergoes regularization through dimensional regularization. In a singular channel, in the case where the interaction is both attractive and sufficiently strong to create a bound state, a pole will emerge below the threshold on the primary Riemann sheet (RS). Conversely, if the attraction is inadequate, the pole will transition to the second RS as a virtual state, persisting below the threshold.

Considering constant contact interactions, saturated by the light vector mesons, and with the coupled-channel effects neglected, in total we obtain a spectrum of 229 hadronic molecules in the charm-anticharm systems. The full results are presented in Ref. [19] and here we discuss on some important systems.

1) The pole positions of the isoscalar $ D\bar D^* $ with positive and negative C-parity are consistent with the molecular explanation of $ X(3872) $^[33] and $ \tilde X(3872) $^[34], respectively. There is a shallow bound state in the isoscalar $ D\bar D $ system, consistent with many theoretical predictions in e.g. Refs. [35−41] as well as the recent lattice QCD result^[42]. Note that the most recent lattice QCD calculations does not find the $ \chi_{c0} $ state corresponding to the $ D\bar D $ bound state^[43–44].

2) The spectrum of the $ \bar D^{(*)} \varSigma_c^{(*)} $ systems is consistent with the molecular explanations of the famous $ P_c $ states^[45−47]: the $ P_c(4312) $ as an isospin-1/2 $ \bar D \varSigma_c $ molecule, and $ P_c(4440) $ and $ P_c(4457) $ as isospin-$ 1/2 $ $ \bar D^* \varSigma_c $ molecules. In addition, there is an isospin-1/2 $ \bar D\varSigma_c^* $ molecule, consistent with the narrow $ P_c(4380) $ advocated in Ref. [24], and three isospin-1/2 $ \bar D^*\varSigma_c^* $ molecules, consistent with the results in the literature.

3) There are two isoscalar $ \bar D^*\varXi_c $ molecules and one $ \bar D\varXi_c $ molecule, which may be related to the recently announced $ P_{cs}(4459) $^[48] and $ P_{cs}(4338) $^[49]. In addition, more negative-parity isoscalar $ P_{cs} $ type molecules are predicted: one in $ \bar D\varXi_c' $, one in $ \bar D \varXi_c^* $, two in $ \bar D^*\varXi_c' $, and three in $ \bar D^*\varXi_c^* $.

4) Instead of associating the $ X(4140) $^[50–51] with a $ D_s^*\bar D_s^* $ molecule like some other works did, our results prefer the $ D_s^*\bar D_s^* $ to form a virtual state. The dip in the invariant mass distribution of $ J/\psi \phi $ measured by LHCb^[52] just at the $ D_s^*\bar D_s^* $ threshold may be produced by this state. A combined fit to the experimental data on the $ 0^{++} $ and $ 2^{++} $ charmonium sector leads to a $ D_s^*\bar D_s^* $ pole as a virtual state^[41]. We also predicted many molecular states in the hidden-charm hidden strangeness systems, most of which are preferred to be virtual states. The peak and dip structures in the invariant mass distribution of $ J/\psi \phi $^[52] may come from these states and a more comprehensive coupled channel analysis of the data is necessary to pin down the signals in such data ¹, as shown in Fig. 2.

Figure 2.  Thresholds of charm-strange meson pairs in the energy range relevant for the $ B^+\to J/\psi \phi K^+ $. Here, $ D_{s0}^* $ denotes $ D_{s0}^*(2317) $, $ D_{s1} $ and $ D_{s1}' $ denote $ D_{s1}(2536) $ and $ D_{s1}(2460) $, respectively, and $ D_{s2} $ denotes $ D_{s2}(2573) $. The data are taken from Ref. [52].

5) The isoscalar $ D^{(*)}\bar D_{1(2)} $ can form negative-parity bound states with both positive and negative C parities, consistent with our previous result^[55]. The $ D\bar D_1 $ bound state is the lowest one in this family, and the $ 1^{--} $ one is consistent with the sizeable $ D\bar D_1 $ molecular component^[56−59] in the $ \psi(4230) $^[60–61]. In addition, we predicted many nearby state with exotic quantum numbers, say, $ J^{PC}=0^{--},\ 1^{-+} $ and $ 3^{-+} $. To check the effects from channel coupling and higher order interactions, we performed a more comprehensive analysis in Ref. [62] with $ D^*\bar D^*\pi $ three-body effects included and it turns out the existence of a $ 0^{--} $ $ D^*\bar D_1 $ molecule is robust.

6) $ \varLambda_c\bar\varLambda_c $ bound states with $ J^{PC}=0^{-+} $ and $ 1^{--} $ are predicted. A fit to the BESIII data on the cross section of $ e^+e^-\to\varLambda_c\bar\varLambda_c $^[63] indicates that the vector $ \varLambda_c\bar\varLambda_c $ bound state should be responsible to the almost flat line shape in the near-threshold region.

7) Light vector meson exchanges either vanish due to the cancellation between $ \rho $ and $ \omega $ or are not allowed in the isovector $ D^{(*)}\bar D^{(*)} $ systems and $ D_s^{(*)} \bar D^{(*)} $ systems. However, the vector charmonium or correlated $ \pi\pi $ exchanges may play important roles^{[30, 64]}, and the $ Z_c(3900,4020) $^[65−69] and $ Z_{cs}(3985)/Z_{cs}(4000) $^{[52, 70]} could well be the $ D^{(*)}\bar D^* $ and $ D_s \bar D^* - D_s^*\bar D$ virtual states. From Fig. 3, we can see that the peak and dip in the $ J/\psi K $ invariant mass distribution from LHCb measurement^[52] coincide exactly with the $ D_s \bar D^* - D_s^*\bar D$ and $ D_s^{*}\bar D^* $ threshold. A coupled-channel analysis on such data may shed light on the nature of the reported $ Z_{cs}(4000) $ and $ Z_{cs}(4220) $ by LHCb^[52].

Figure 3.  Thresholds of $ D_s\bar D^* $ and $ D_s^* \bar D^* $ in the energy range relevant for the $ B^+\to J/\psi \phi K^+ $. The data are taken from Ref. [52].

When including exchange of light vector mesons, our results are generally consistent with those obtained from more complete treatments of the one-boson exchange model by solving the Schrödinger equation.

Examples include the binding energies of isoscalar $ D\bar D $^[35] and $ D\bar D_1 $^[55] bound states from ρ and ω exchanges, and a similar pattern of molecular states related to the $ X(3872) $ in Ref. [71], where light vector exchange was found necessary to bind $ D\bar D $ together. $ \bar D \varSigma_c $ bound states corresponding to the $ P_c(4440) $ and $ P_c(4457) $ were obtained through one-boson exchange in Ref. [72], and the degeneracy of the two states with $ J=1/2 $ and 3/2 is lifted by pion exchange and higher order contributions from ρ and ω exchanges. It is worth noting that some systems may have contact terms receiving important contributions from scalar meson exchanges. Among the 229 predicted structures, only a few have been studied in detail due to energy limitations in current experiments. To find more states in the predicted spectrum, we need higher statistics and data from other experiments, such as prompt production at hadron colliders, PANDA, electron-ion collisions, and $ e^+e^- $ collisions above 5 GeV at super tau-charm facilities.

Considering the constant contact interactions saturated by the light vector meson exchange while disregarding the effects of coupled channels, we obtain a spectrum of 124 hadronic molecules in total in the charm-charm systems. At least the same number of molecules are expected to exist for each of the charm-bottom and bottom-bottom systems, given that forming a bound state is more feasible with similar attraction strengths due to the heavier reduced masses; Furthermore, there may be additional molecules beyond this estimate, especially if the ground state is significantly bound, indicating the possible existence of excited states. The full results are presented in Ref. [19] and here we discuss on some important systems.

1) The isoscalar $ D^*D^{(*)} $ systems exhibit an attractive interaction due to light vector meson exchange, leading to a total potential that places these systems on the verge of forming near-threshold molecules. By applying a reasonable cutoff to regularize the loop integral, the binding energy of the $ I(J^P)=0(1^+) $ $ DD^* $ system aligns with the observation of the double charm tetraquark $ T_{cc}^+ $ by LHCb^{[73, 74]}. This configuration is commonly interpreted as the $ I(J^P)=0(1^+) $ $ DD^* $ molecule^{[75, 76]} and in such case, more similar states with $ I=0 $, such as $ D^*D^* $, $ D^{(*)}D_{1,2} $, and $ D_{1,2}D_{1,2} $ are predicted to exist.

2) Note that the attraction of the $ I=1/2 $ $ D^{(*)}\varSigma_c^{(*)} $ and $ D_{1,2}\varSigma_c^{(*)} $ systems arising from the exchange of light vector mesons is more potent for those in double-charm $ D^{(*)}\varSigma_c^{(*)} $ and $ D_{1,2}\varSigma_c^{(*)} $ systems. It is reasonable to anticipate the presence of such double-charm hadronic molecular states by considering the well-known $ P_c $ states being essentially hadronic molecules composed of $ \bar D^{(*)}\varSigma_c^{(*)} $. A similar argument can be extended to the channels involving $ D^{(*)}\varXi_c^{(\prime*)} $ and $ D_{1,2}\varXi_c^{(\prime*)} $, particularly in light of potential experimental confirmation of $ P_{cs} $ states.

3) Within our simple model, we anticipate that only the isoscalar $ \varSigma_c^{(*)}\varSigma_c^{(*)} $ systems in the double-charm di-baryon sector could potentially form bound di-baryon states.

By examining the shape of the invariant mass distribution in a lower channel close to the threshold of a coupled higher channel, we expect to identify nontrivial structures in the near-threshold region for all heavy-antiheavy hadron pairs that exhibit attractive interactions at the threshold. These structures might appear precisely at the threshold as threshold cusps if the attraction is not strong enough to create a bound state, or they could manifest below the threshold in the presence of a bound state. It is important to notice that these structures may not always appear as sharp peaks but could also be observed as dips in the invariant mass distributions, depending on how large are the production rated of each channel. These characteristics align with various experimental findings, suggesting the likelihood of discovering more near-threshold structures in experimental studies.

To proceed, we systematically calculated the spectrum of hadronic molecules comprising a pair of charmed and anticharmed hadrons as well as a pair of charmed hadrons. This was accomplished through the utilization of S-wave constant contact potentials that are saturated by the exchange of vector mesons. Our study encompasses all S-wave singly-charmed mesons and baryons, as well as the $ s_\ell=3/2 $ P-wave charmed mesons, by incorporating HQSS, chiral symmetry, and SU(3) flavor symmetry in establishing the coupling between charmed heavy hadrons and light vector mesons. It is important to view our predicted spectrum as a primary approximation of the heavy-antiheavy and heavy-heavy molecular spectrum, which outlines general patterns of such hadronic molecules. The approximation we have used may introduce significant quantitative adjustments in the numerical results. Specifically, the following limitations should be noted:

1) The momentum-dependent terms, encompassing both spin-dependent and spin-independent contributions, have not been taken into account.

2) The impact of coupled channels has been overlooked in this work. In certain scenarios, these effects could play a crucial role in the development of states close to the energy threshold.

3) Double-charmed meson-baryon molecules have the ability to interact with regular double-charm baryons as well as channels that consist of a double-charm baryon and a light meson.

4) The interactions discussed here are limited to the leading order in the $ 1/N_c $ expansion with $ N_c $ the number of colors. Consequently, we have overlooked the Okubo-Zweig-Iizuka suppressed interactions that could potentially play a role in resolving the degeneracy within the system characterized by varying total spins. Furthermore, we may also miss some states that arise from these OZI suppressed interactions, such as the potential double-charmonium states discovered in experiments^[77−79] and the isospin partner of $ X(3872) $ predicted recently in Ref. [80].

Hence, while we anticipate that the spectrum outlined here will illustrate a general pattern of the hadronic molecules formed of a pair of charmed mesons and/or baryons, individual systems might deviate quantitatively from the prediction.

The potentials in the bottom sector are the same as those in the charm sector, assuming nonrelativistic field normalization, owing to heavy quark flavor symmetry. It is anticipated that a similar number of molecular states will be present in analogous systems. Given the higher reduced masses of hidden-bottom systems, virtual states in the charm sector will approach the thresholds or potentially become bound states in the bottom sector, leading to more profound binding of bound states in the charm sector within the bottom sector.

Acknowledgments I am grateful to Feng-Kun Guo and Bing-Song Zou for their insightful discussions and fruitful collaborations on these projects.