Coupled-channel description of charmed heavy hadronic molecules within the meson-exchange model and its implication

Motivated by the first observation of the double-charm tetraquark $T_{cc}^+(3875)$ by the LHCb Collaboration, we investigate the nature of $T_{cc}^+$ as an isoscalar $DD^*$ hadronic molecule in a meson-exchange potential model incorporated by the coupled-channel effects and three-body unitarity. The $D^0D^0\pi^+$ invariant mass spectrum can be well-described and the $T_{cc}^+$ pole structure can be precisely extracted. Under the hypothesis that the interactions between the heavy flavor hadrons can be saturated by the light meson-exchange potentials, the near-threshold dynamics of $T_{cc}^+$ can shed light on the binding of its heavy-quark spin symmetry (HQSS) partner $D^*D^*$ ($I=0$) and on the nature of other heavy hadronic molecule candidates such as $X(3872)$ and $Z_c(3900)$ in the charmed-anticharmed systems. The latter states can be related to $T_{cc}^+$ in the meson-exchange potential model with limited assumptions based on the SU(3) flavor symmetry relations. The combined analysis, on the one hand, indicates the HQSS breaking effects among those HQSS partners, and on the other hand, highlights the role played by the short and long-distance dynamics for the near threshold $D^{(*)}D^{(*)}$ and $D^{(*)}\bar{D}^{(*)}+c.c.$ systems.


I. INTRODUCTION
One of the critical issues about the non-perturbative property of quantum chromodynamics (QCD) is to what extent it allows the existence of the so-called "exotic hadrons" of which the constituent contents are beyond the conventional quark model (i.e.mesons are composed of q q and baryons of qqq [1][2][3][4][5]).Such exotic hadrons include glueballs, hybrids, multiquarks, and hadronic molecules, etc, and their existences may serve as a unique probe for understanding the non-perturbative property of QCD.Since the discovery of X(3872) by the Belle Collaboration in 2003 [6], there have been a large number of exotic candidates observed in experiments (see e.g.Refs.[7][8][9][10][11][12][13][14][15][16] for recent reviews).Interestingly, most of these observed states are heavy flavor states and located in the vicinity of some relative S-wave thresholds.It makes the hadronic molecule picture a natural solution for their nature, and also allows the implementation of effective field theory (EFT) approaches in the description of the near-threshold dynamics.Such a phenomenon is very similar to that for the deuteron system, where the EFT approach originally proposed by Weinberg [17] have been greatly developed during the past decades [18][19][20][21].
From the perspective of the EFTs, the non-perturbative property of QCD can be learned from different aspects at different low energy scales.In the light-quark sector, the light meson-meson interactions can be described by the chiral perturbation theory (ChPT), which is constructed by the pseudo-Nambu-Goldstone boson fields that emerge from the spontaneous symmetry breaking of chiral symmetry.While in the heavy-heavy sector, the interactions between two heavy hadrons are well constrained by the heavy quark symmetry (HQS) and heavy quark spin symmetry (HQSS).Due to the large heavy quark masses (m Q ≫ Λ QCD ) the mid and short-range interactions can be described by only several low energy constants (LECs) for different isospin channels.However, these LECs still cannot be well-determined due to the lack of experimental data and taking into account also the uncertainty O(Λ QCD /m Q ) if m Q is not large enough.
From another view, the hadron-hadron interactions can also be described by the potentials generated by the meson exchanges, such as scalar (σ, κ), vector (ρ, ω, J/ψ), pseudoscalar mesons (π, η, η c ), etc.They are similar to the meson exchange potentials between nucleons inside a nucleus and can mimic the underlying quark-exchange processes.With the decay constants determined by the on-shell heavy-meson decays, these LECs are hypothesized to be saturated by matching the ChPT operators expansion with those meson-exchange potentials [22,23].
With the interactions mentioned above, the dynamical mechanism between hadron-hadron pairs can then be investigated and the complete scattering amplitudes can be accessed by making a proper analytical continuation of the a E-mail address: qiulin@ihep.ac.cn b E-mail address: gongchang@ccnu.edu.cnc E-mail address: zhaoq@ihep.ac.cn on-shell ones into the complex plane with the exception of finite singularities determined by the particle-exchange kinematics [24,25].Then the resonance peaks observed in the invariant mass spectrum may be derived from the poles of the scattering amplitudes on the complex plane, which has been proved successful in the description of many low-energy meson-meson scatterings [26,27] and meson (nucleon)-nucleon scatterings [28,29] in the light-quark sector.As to the very long-range potential mediated by the light pion exchange, it brings into the breaking of HQS/HQSS, and due to the on-shell decay D * → Dπ the two-hadron composite systems are usually quasi-bound/virtual-states-like resonances.Strictly, the multi-hadrons involved processes can be learned by solving the Faddeev-type equations [30].But in practice they are often reduced to an effective two-body Lippmann-Schwinger equation with the assumption that the two-body interaction proceeds via an isobar [31,32].In such a way, the three-body unitarity is reserved and the poles of these hadronic molecules acquire a certain imaginary part which arises from these sub-threshold cuts.
In 2021, the LHCb Collaboration announced the first observation of a double-charm tetraquark T + cc (3875) (ccū d) in the D 0 D 0 π + invariant mass distribution [33] and its Breit-Wigner parameters are, δm BW = −273 ± 61 keV/c 2 , Γ BW = 410 ± 165 keV, with δm BW the mass relative to the nominal D 0 D * + threshold and Γ BW the width.Later, a unitarized analysis by LHCb Collaboration [34] suggests it is a pole on the second Riemann sheet, where δm pole also refers to the real part of its pole relative to the D 0 D * + threshold and Γ pole is twice of the absolute value of the imaginary part of its pole.The partial amplitude analysis and the absence of its isospin partners in D + D 0 π + channel imply its quantum numbers as IJ P = 01 + .Besides, the fact that approximately 90% of the D 0 D 0 π + events contain a genuine D * + meson reveals its main components as D 0 D * + , and it should couple strongly to D + D * 0 and DDπ due to the proximity among those channel thresholds.Early theoretical studies of the double heavy-flavor exotic states can be found in the literature [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49].
Different from the double-charm system, the charmed-anticharmed systems have already been studied broadly in both experiment and theory.Although the nature of these resonances is still heatedly debated due to the more complicated analytical structures.In our work, we first study the near-threshold dynamical nature of T + cc as an isoscalar DD * hadronic molecule within the meson-exchange model including both coupled-channel effects and the three-body unitarity.Then, we try to explore other possibilities of hadronic molecules composed by charm-anticharm heavy mesons such as X(3872) and Z c (3900) which can be regarded as partners of T + cc with respect to flavour symmetry and charge conjugate transformation beyond the HQS/HQSS in a more general basis.In this sense, our combined analysis of these systems tends to, on the one hand, demonstrate the relations based on the HQS/HQSS, and on the other hand, manifest the effects caused by the HQS/HQSS breaking in the D ( * ) D( * ) and D ( * ) D ( * ) interactions.
This paper is organized as follows: In Sec.II, we discuss the potentials constructed by the one-boson-exchange (OBE) model and coupled-channel formalism, and the eigen wavefunctions satisfying the Hamiltonian.In Sec.III, we present our numerical results and discussions.A brief summary is given in Sec.IV.

II. FRAMEWORK A. One-boson-exchange model and potentials
The OBE model describes the hadron interactions similar to that for the nuclear force.The exchanged bosons, which have different quantum numbers and masses, account for a lot of features of the nuclear elements.While the boson exchanges are actually the quark exchanges between hadrons, the OBE model can be regarded as a leading order approximation of the strong QCD in the non-perturbative regime.Within the OBE model the dynamics for the interacting hadrons are generally fine-tuned by a physical cutoff Λ (0.3 ∼ 1.2 GeV) by matching the theoretical calculations to the experimental data.The cutoff indicates a typical scale beyond which the boson exchange scenario does not hold anymore.The relevance of these dynamical ingredients can be learned by the so-called Weinberg organizational principle [18][19][20][21] which efficiently collects the soft/hard scales and counts the operators order by order.
In the heavy meson-(anti) meson sector [70], it can be learned that the one-pion-exchange (OPE) and coupledchannel effects (i.e. the DD − DD * − D * D * couplings for the double-charm system), are both next-to-leading-order (NLO).In contrast, the leading-order (LO) contributions are from the light vector meson exchanges, which are treated as dynamical gauge bosons of hidden local symmetries [71][72][73][74][75].At this moment, we limit our analysis in the sector of non-strange charmed meson systems.We note that in the bottom sector the large masses of B ( * ) mesons can bring in a large typical relative momentum for the B ( * ) B (B ( * ) B) systems at the mass region of the B * B * (B * B * ) threshold though the threshold differences among B ( * ) B ( * ) or (B ( * ) B( * ) ) channels are rather small, i.e. δ < m π with δ ≡ M B * − M B ≃ 50 MeV.As the result, the OPE contribution is no longer perturbative and can be promoted to the LO [76].This will slow down the convergence of the EFT.For composite systems with strangeness, the SU(3) flavour symmetry becomes the approximate one and the kaon-induced interactions are very subtle due to lightness of kaon mesons.
For composite system DD * in an S-wave, the wave functions with quantum numbers J P = 1 + with I = 0, 1 can be constructed as, It should be noted that the systems of D 0 D * 0 (I = 1, I 3 = −1) and D + D * + (I = 1, I 3 = +1) are just partners of the DD * isospin triplet with (I = 1, I 3 = 0).In view of the hidden local symmetries the isoscalar system is found to be attractive, and the isovector one turns out to be repulsive with the light vector meson exchanges.Thus, we consider only the above wave functions for DD * and charge-neutral ones for the charmed-anticharmed system.Note that the full eigenstate must include all the degrees of freedom.By solving the two-body problem, one obtains the eigenstate for this two-body system, but would not be able to distinguish D * and D. The eigenstate describes the relative motions of the two-body system quantum-mechanically in a relative spatial separation of r (corresponding to the momentum transfer between the two constituents in the momentum space) without knowing which one is which.Similar considerations were also discussed in Ref. [84].In our recent work [85], we have made a detailed survey of the wave functions and potentials of the heavy hadronic molecules in accordance with their quantum numbers and we briefly summarize the basic ingredients here for reference.This issue is also noticed by Ref. [57].
Moreover, we emphasize that the total wave function, which is symmetric under space, isospin and spin group O(3)⊗ SU I (2) ⊗ SU S (2) for bosons, is crucial for properly introducing the meson-exchange potentials, especially the pionexchange ones from the u channel.Since the full pionful interaction requires the inclusion of loop corrections [86,87] and the pion-exchange potential would not account for the short-distance dynamics with r << 1/m π , different parametrizations of the short-distance potential can result in different behaviors of the nonperturbative OPE to be either attractive or repulsive [70].This also calls for a unified approaches for the D ( * ) D( * ) and D ( * ) D ( * ) systems.Otherwise, the conclusion about the importance of OPE could become self-contradicting and even misleading.
For the S-wave D D * composite system, the charge-neutral wave functions for different isospins and charge conjugations can be constructed in a similar way, Since for many of those hadronic molecules, the binding energy is far smaller than the threshold difference between coupled channels and thus the isospin breaking effect should be taken into account.Therefore, we define two channels in the particle basis as in Ref. [81]: The physical state [81,88,89]) taking into account the coupled-channel effects.Similarly, we define two channels for D D * : The near-threshold potentials can then be evaluated by the static approximation [90]: with M the scattering amplitude of the process ab → cd and m i/f the masses of initial or final state mesons.For convenience, we first define the following functions: where ϵ i is the initial polarization vector of D * ( D * ), ϵ * f is the final one, m ex is the mass of the exchanged meson, and q 0 /q are the zero-th/three-vector components of the transferred momentum q µ = p ′µ − p µ with p ′µ /p µ the momenta of the final/initial mesons.Interestingly, q 0 could be larger than the mass of the exchanged pion for some processes like It will lead to a logarithmic divergence in the projected partial waves and should be properly treated by a full consideration of the three-body unitarity [91,92].
The detailed potentials for DD * → DD * (as depicted in Fig. 1) are listed as follows: Similarly for D D * , To saturate the LECs with the meson-exchange potentials, one matches the contact-range operator expansion with those potentials generated by the light meson exchange [93].We illustrate this pattern with OPE by decomposing it into two terms, where with q the unit of q.Because of the perturbative contribution by the pion-exchange diagrams, the tensor term is negligible and we only consider the S-wave component in our work.For the central term, the first Dirac delta term actually is unphysical since the dynamics of δ(r) is far beyond the reach of the pion scale 1/m π .It is generally parametrized by a finite size ∼ 1/Λ with Λ > m π .In our approach, we parametrize this term along with the heavier pseudoscalar (P ) and vector (V ) meson (m P,V ∼ Λ) exchange potentials via the scale-dependent counter terms C P,V , which are not fully captured by the OBE model [94].Thus, we adopt the following substitution rules [95], with µ 2 i = m 2 i − q 0 2 and 0 ≤ C P/V ≤ 1.For very heavy exchanged mesons (e.g.J/ψ, η c ), we just ignore the spin-dependent terms Y ex , Z ex as they vanish at the threshold (since the Dirac delta terms are fully reserved), but keep the resumed Yukawa terms X ex .We collect the wave functions and potentials of other heavy double-charm or charmed-anticharmed systems in Appendix A.

B. Lippmann-Schwinger equation and three-body interactions
The near-threshold coupled-channel dynamics can be studied by solving the nonrelativistical Lippmann-Schwinger equation (LSE), where E is the energy in the initial c.m. frame, V αβ (p, k; E) is the potential from the β-th channel to the α-th channel.
To regularize the ultra-violet (UV) divergence we introduce a hard cutoff into the LSE by replacing the potentials with, with Θ the Heaviside step function and Λ the cutoff parameter.We refrain from the use of any Gaussian or monopole form factors because of their complexities in the undermentioned analytical continuation [96].G δδ (q; E) is the twobody propagator of the δ-th channel which can be written as (in the nonrelativistical limit), with m i the mass of the constituent particle, µ δ the reduced mass of the δ-th channel and Γ δ (q; E) the width contribution of the δ-th channel.As we only focus on the near-threshold phenomena (e.g.3.873 ∼ 3.877 GeV for the D 0 D 0 π + mass spectrum), the above nonrelativistic approximation should be reasonable.Since the width of D * is comparable with the width of the corresponding hadronic molecules, namely T + cc and X(3872), the width part of the propagator Γ δ will have a significant effect on the pole position.In most of the literature, the constant width treatment for D * might be appropriate as argued in Ref. [97].But the width of the composite state could be overestimated by two times, compared with the one in which the full three-body unitarity is considered [91].In the case of a constant width for D * , the Riemann sheets can be classified according to the sign of the imaginary parts of the on-shell momenta (taking two channels D 0 D * + and D + D * 0 as example): with µ 1 , µ 2 the reduced masses of the D 0 D * + , D + D * 0 channels, respectively.Namely, RS-I : Apart from the LO tree diagrams at the order Q −1 , the NLO diagrams involving the pion exchange which can contribute to the decay width of composite particle DD * ( D * ), is depicted in Fig. 2. They are suppressed by one power of Q as analyzed with an effective Lagrangian derived from heavy-hadron chiral perturbation theory (HHχPT) [98].While the NLO contribution is also found to be dominant by the contact potential, the wave function receives the renormalization, especially the imaginary part, by evaluating the above cut diagrams which are characterized by the Dπ loop and one-pion-exchange.Regardless of the full-pion calculation, the nonperturbative nature of the OPE was justified by the opposing argument about the calculation of one-loop diagrams using the dimensional regularization scheme [86] and the sharp cutoff scheme [87] which leads to the conclusion that the OPE potential in an EFT is only well-defined in connection with a contact term.Moreover, the width of X(3872) with the effects of three-body cuts and its quark-mass-dependence were also validated in the EFT treatment with both perturbative pions [98,99] and nonperturbative pions [91].In summary, the DD * scattering in three-body unitarity must include the full D * propagators and the pion-exchange diagrams in the on-shell renormalization scheme.
In our approach since the dominant contact part in the OPE is completely removed, the perturbative nature of OPE is retained if the physical cutoff and other counter terms are fine-tuned to the full line shape.We leave the discussion about its analytical treatment in Appendix B. On the other hand, the width of D * should be energy-dependent due to the self-energy from Dπ loop as done in Refs.[52,91,92,95].For D * → Dπ, the width contribution to the complex mass of D * is, where In the moving frame of D * (q), the energy available for D * gets reduced since the denominator of the Dπ propagator is now, where ω i = m 2 i + l 2 is the energy of the constituent particle in the static c.m. frame of D * (q = 0) with l the relative three-momentum between D and π in the Dπ loop.In the framework of Time-Ordered-Perturbation-Theory (TOPT), the static energy of D * in the initial c.m. frame can be approximated by with m D * the bare mass of D * .We denote the p cm solved from the above process by: p cm = F(E, q, m 1 , m 2 , m 3 ) and then the energy-dependent width of D * + /D * 0 is (analogous to Refs.[91,92]), For the purpose of searching for poles, the function F(E, q, m 1 , m 2 , m 3 ) has to be analytically and properly continued by the technique of redirecting the branch cuts and contour deformation.We only present the result in this section and refer to Refs.[92,96,100,101] or Appendix B for details, C. Line shape analysis of the D 0 D 0 π + mass spectrum We solve the Lippmann-Schwinger equation in the particle basis, and the LO Lagrangian to describe the interaction between T + cc and D/D * is, where symbol ∓ corresponds to the isoscalar/isovector DD * system, and g XDD * is the coupling constant and approximated by a constant due to the narrow D 0 D 0 π + invariant mass range ∼ 4 MeV.Note that the contribution from the D wave is neglected due to the suppression of the centrifugal barrier.Thus, the transition T + cc → D 0 D * + can be evaluated by (as depicted by Fig. 3), with p the three momentum of D 0 in the c.m. frame of T + cc .The coupling constant g XDD * can be absorbed into the overall factor N and thus is set to unit.The decay width of T + cc → D 0 D 0 π + is (see Fig. 4), where the two terms come from the symmetry of two D 0 mesons in the final state with q π (q π ) the magnitude of three-momentum of pion in the c.m. frame of D * + [π + (p 2 )D 0 (p 3 )] and p (p) the magnitude of the D 0 three-momentum produced at the first vertex in the initial c.m. frame, namely, where λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz is the Källén function, s 12 = (p 1 + p 2 ) 2 , s 23 = (p 2 + p 3 ) 2 are the invariant masses squared, and the integral limits are determined by the Dalitz boundary, with m 1 , m 2 , and m 3 the masses of D 0 , π + , and D 0 , respectively, and and . In order to compare with the experimental data, the decay width function above should be convoluted with the mass resolution function of the LHCb detector, which are modulated by the sum of two Gaussian functions [54], where ∆E = 200 keV is the bin width; β 1 = 0.778, β 2 = 0.222, σ 1 = 1.05 × 263 keV, and σ 2 = 2.413 × σ 1 are taken from the LHCb analysis [33,34].

III. RESULTS AND DISCUSSIONS
A. Line shape of D 0 D 0 π + mass spectrum Proceeding to the numerical analysis, we list the parameters to be fitted by the line shape data: the physical cutoff Λ, the counter terms C P/V , and the trivial overall factor N .It is obvious that Λ and C P/V are correlated but we find that they can be strictly constrained by fitting the line shape of the D 0 D 0 π + mass spectrum.Moreover, C V has a more sensitive impact on the pole position and line shape while the dependence of C P is relatively moderate due to their coupling difference and the energy scale difference between the pseudoscalar and vector mesons in such systems.
In order to investigate the role played by the OPE mechanism and the width effects of D * , we consider three fitting schemes to fit the D 0 D 0 π + mass spectrum with Eq. ( 37), i.e.
3. Scheme III: OBE potentials with an energy-dependent D * width which incorporates with the three-body unitarity (i.e. the OPE is also properly considered).
In the above three schemes we consider a baseline fit by fixing C P/V = 0 to compare with the fitting results with C V fixed with well-chosen values from 0.0 to 1.0.
In Table .I the fitted parameters (when C V = 0.0) and the extracted pole positions are presented.The fitted line shapes in the corresponding schemes are presented in Fig. 5. Roughly speaking, these fitting results are comparable and reasonably good.But by detailed comparisons we can still learn some crucial information concerning the underlying dynamics: • Since the contact-range term of the OPE has been parametrized out by parameter C P which is then set as C P = 0, it means that the OPE will only account for the long-distance behavior of the wavefunction, while the short-distance behavior will be parametrized by C V .As the consequence, the cutoff Λ in the three schemes is consistently in the same range.We mention that if C P ̸ = 0 the fitted values for Λ can be very different for these three schemes.This can be compared with similar findings in Ref. [92].
• Although the OPE contribution is subleading, it still plays an important role in fitting the line shape and extracting the pole position and width, which can be seen by comparing between Scheme-II and Scheme-III.Moreover, the width at the bound state pole including three-body unitarity (Scheme-III) is consistent with the unitarized analysis by LHCb [34].
• Also, it shows that the implementation of the energy-dependent D * width will reduce the pole width by a half.This indicates the importance of the three-body unitarity to some extent.
As mentioned earlier, there exist strong correlations between the counter terms C P/V and Λ.Thus, we make a survey of their correlation by looking at the variations of Λ and χ 2 with C P being fixed to zero and C V taking values within [0, 1].The correlation is demonstrated in Fig. 6.It shows that, as C V increases, the Λ increases with different increasing rates in the small (physical) and large (unphysical) momentum regions (see Appendix.C for a detailed explanation).The cutoffs in the three schemes, especially in the physical region (≤ 0.5 GeV), are consistent with each other due to the NLO roles played by the three-body cuts.FIG.5: Baseline fitting results of line shape of D 0 D 0 π + , convoluted with mass resolution function.-I) in Scheme-I is defined by the two-body branch point given by Eq. ( 23) and the second Riemann sheet (RS-II) in Scheme II-III is the most important unphysical Riemann sheet accessed by Eq. ( 29).The uncertainty of Λ is evaluated by χ 2 fitting and propagates to pole positions.FIG.6: (Left) Dependence of Λ on C V in these three schemes.(Right) The fluctuations of χ 2 fitting on C V in these three schemes.The slight instability of χ 2 in Scheme-III is caused by the mild dependence of the fitting scheme on parameter ω used in the contour deformation.But one can read from the left panel that such a fluctuation does not cause significant deviations of Λ from its smooth correlation with C V .
For isoscalar systems excluding 1 +− D D * and 2 ++ D * D * (see Appendix A), the larger C V is, the larger repulsive contact potential will be introduced and thus the larger Λ is required to retain more short-distance potentials.Meanwhile, for a relatively large Λ, more contributions from higher-order terms, such as O(q 2 /(q 2 + µ 2 π )) of OPE, may manifest themselves as shown from the small but sizable deviation (due to the weak scaling of pion-exchange) between Scheme-II and Scheme-III.Nevertheless, the reasonable range of the cutoff parameter Λ is set empirically below 1.2 GeV.
Actually, the stability of the χ 2 for the whole range of C V with the corresponding range of Λ in Fig. 6 indicates the model-independent feature in this analysis after the proper treatment of the short and long-distance dynamics and the three-body unitarity [102].This also reflects the matching between the OBE and the EFT approach (see e.g.Ref. [92]).Such a feature can be further investigated in the study of the D * D * (I = 0) system as the HQSS partners of T + cc .On the one hand, we expect that the HQSS relation should provide some constraints on the parameters fixed by the T + cc line shape.On the other hand, we also anticipate deviations which are the manifestations of the HQSS breaking effects.

B. Possible hadronic molecules in the nonstrange charmed sector
For the double-charm system, the only possible hadronic molecule except T + cc is the isoscalar 1 + D * D * and the contact potentials between them are strictly correlated by the heavy-quark spin symmetry.Namely, in the heavy quark limit, we have Thus, the parameters learned from T + cc can make a strong constraint on the binding energy E B of the isoscalar D * D * .We plot the dependence of E B on the aforementioned C V (Λ) in the three schemes as in Fig. 7 and the E B is found to be about 2 ∼ 9 MeV for reasonable parameters (C V < 0.7, Λ < 1.0 GeV), which is consistent with the result in Ref. [53] (but the isoscalar D * D * is only found to be a bound state here).As to the charmed-anticharmed system, the nonperturbative short-distance interactions can be very different from the double-charm ones, such as the difference between the quark annihilation of cc and quark rearrangement of cc.Such dynamical differences cannot be fully compensated by the meson-exchange model, and neither by an EFT approach given the HQSS breaking is unavoidable.To be consistent with the experimental results, it is found that C V tends to be much larger under the same Λ when C P is still fixed as zero.For example, the binding energy of D * D * (I = 0) would be found around 100 MeV if the parameters from T + cc adopted.Such a value has obviously overestimated the binding since the binding momentum |p| would amount to hundreds of MeV when its binding energy is just 10 MeV.Meanwhile, it will be questionable for the Born approximation since it would not hold for such deeply bound states.On the other hand, taking into account that the isoscalar D D * with C = + is found to form a bound state for a large C V , but only a virtual state for a small C V , the existence of X(3872) allows us to limit our discussions within 0.35 < Λ < 0.65 GeV and a relatively large 0.7 < C V < 1.0 below though we present the results within a wide range of C V in Fig. 7 to show the Λ dependence [103].The LSE solutions for the isoscalar D ( * ) D( * ) are listed in Tab.II.
TABLE II: Pole positions (in MeV) of the isoscalar D ( * ) D( * ) .The notation (rs, E B ) represents the binding energy E B on the first/second Riemann sheet.The E B boundary for each Λ is evaluated with C V = 0.7, 1.0, respectively.See the context for further explanations.
Since there is no spin-dependent terms for D D, the above Λ boundary is chosen partly by comparison with another calculation by dimensional regularization (DR) in Ref. [53].While for D * D * , there is a mass splitting for different J due to the inclusion of spin-dependent terms and we note that the unacceptably large numerical binding energy 117.7 MeV for the 2 ++ D * D * system is calculated with the unphysical value of C V = 1.0.In such a case the full reservation of the Dirac Delta term will introduce an extremely large attractive potential as mentioned in the last subsection.Though more profound considerations, such as higher partial waves, coupled-channel effect, and four-body unitarity, etc, are needed, the above result shows the signal for the existence of the isoscalar hadronic molecules composed of D * D * , especially the 2 ++ tensor state.Despite the fact that there is no lower channels than the D D threshold in strong decays, the existence of an isoscalar D D has been widely predicted by phenomenological studies [104,105] and by LQCD calculation [106].Besides, the role of the isoscalar D D state can be significant in the study of its HQSS partners like X(3872) by the loop enhancement if there exists a near-threshold bound/virtual state [107,108].As for the isovector D ( * ) D( * ) , we found no bound/virtual states here but the observed hadronic molecule candidate Z c (4020) ± [109,110] is actually generated by more complicated mechanisms, similar to Z c (3900) to be discussed later.
For the D D * system, the three-body unitarity has to be taken into account.Interestingly, it seems that there should exist a counterpart of X(3872) with negative C-parity, denoted as X(3872), which has been studied in the literature [111][112][113][114].In Fig. 8, we plot the pole positions of these two states with respect to Λ and C V .Similarly, the binding energy of X( 3872) is estimated to be from hundreds of keV to a few MeV along with an imaginary part comparable with that of T + cc .One notices the different dependence behaviors of the width and binding energy E B on the short-distance parameter C V for X(3872).Namely, with the increase of C V the value of |E B | drops and the width (defined as twice of the absolute value of the imaginary part of the pole position) increases.In particular, with C V → 1, the width increases up to 40 ∼ 60 keV where E B is only hundreds of keV.This behavior is consistent with the results of Refs.[91,98].In the upper panels of Fig. 8 one also sees the correlation of the width and binding energy with the cutoff parameter Λ with the fixed C V .With the smaller value for Λ, one obtains a relatively shallower binding and relatively larger width.Note that with the fixed C V the smaller value of Λ corresponds to the smaller contributions from the short-distance dynamics.The driving mechanism for the binding of X(3872) due to the short-distance dynamics is actually highlighted.
In contrast with X(3872), the C = −1 isoscalar X(3872) is found to be more bound with a slightly smaller imaginary part.Actually, if one only looks at the solid lines for both X(3872) and X(3872) in Fig. 8, their values are quite close to each other.However, with the increase of the Λ value the binding energy and total width of X(3872) both increase.As mentioned earlier that the larger value of Λ corresponds to larger contributions from the short-distance dynamics, it suggests that the coupling strength for X(3872) to D D * + c.c. increases.Such a change, on the one hand, leads to a larger binding energy, and on the other hand, overtakes the decrease of the phase space to produce a larger decay width.With Λ ∼ 0.4 GeV favoured by X(3872) (by comparing the extracted pole mass with the value from the Particle Data Group (PDG) [115]), the binding energy of X(3872) is about 10 MeV, which is very close to the result M X(3872) = 3860.4± 10.0 MeV/c 2 reported by the COMPASS Collaboration [111].However, it should be cautioned that for such a deeply bound state, the formula of the D * width must be revised regarding the small mass difference between D * and Dπ, and the large binding energy up to 10 MeV.To be brief, comparing the behaviors of X(3872) and X(3872) with each other one sees the crucial role played by the short-distance dynamics.Meanwhile, one sees the apparent deviations from the HQSS.
For the isovector D D * system, the interactions from the ρ and ω meson exchanges actually cancel with each other.Without additional mechanisms this system cannot form bound state through the residual meson-exchange potentials (i.e.σ, π, η, J/ψ, etc), and it turns out more likely to be a virtual state or resonance (see e.g.Refs.[116][117][118][119][120] and the review of Ref. [7]).This can be regarded as being consistent with the experimental observations of Z c (3900) [121][122][123].In particular, the observation of Z c (3900) seems to be strongly correlated with the production of Y (4260) for which the presence of the so-called "triangle singularity" can provide a natural explanation for its enhanced production rate [124,125].

IV. SUMMARY
Motivated by the observation of the double-charm tetraquark molecule candidate T + cc by the LHCb Collaboration, we present a coupled-channel analysis of the D 0 D 0 π + mass spectrum within the meson-exchange model including isospin breaking effect.It is found that the dynamical details, and the two-body or three-body unitarity, have indeed played an important role in the description of the underlying dynamics for such composite states with unstable constituents.The well-described line shape allows us to extract the pole information about T + cc at a rather high precision and it confirms the nature of T + cc as an isoscalar DD * hadronic molecule which can be regarded as a milestone in the study of hadronic molecules and hadron spectroscopy.Moreover, by implementing the HQSS it provides a stringent constraint on some of those crucial dynamical aspects for the D ( * ) D( * ) systems which allows us to compare with the EFT results and examine the HQSS breaking effects through the combined analysis.In particular, the existence of another double-charm hadronic molecule candidate D * D * (I = 0) is predicted in the same framework.
In the combined analysis the interactions between the charmed mesons are hypothesized to be saturated by the meson-exchange potentials.This allows us to relate the OBE potentials for the double-charm system with the charmed-anticharmed system.Although it should be recognized that the short-distance dynamics between the charm quark rearrangement and charm-anticharm annihilations are actually different, we expect that these two systems will still share some crucial aspects of the dynamics via the OBE potentials.It allows us to establish the relations between the double-charm and charmed-anticharmed systems, and extract the binding conditions for these systems, especially for the isoscalar ones.For isovector states such as Z c (3900) and Z c (4020), we find that the light vector mesons potentials via the exchange of ρ 0 and ω actually cancel each other.The residual potentials turn out to be insufficient for binding.It suggests that Z c (3900) and Z c (4020) should behave more likely to be threshold-enhanced virtual states, resonances, or even quasi-bound states depending on additional mechanisms introduced.Besides, the correlation of the Z c (3900) and Z c (4020) productions with the triangle singularity mechanism make it necessary to include the production mechanisms in the detailed analysis.While these are still non-trivial issues to be addressed, we leave them to be studied in future works.
+ ] is static, i.e., q = 0.The three-body cut is redirected from the RHC on the real axis to the direction parallel to the negative imaginary axis.
Meanwhile, the integral path of q shall be deformed and we adopt the parametrization used in Ref. [96], namely, With ω ∼ 0.1 GeV and V 0 ∼ −0.1 GeV, the "spectator momentum contour" (SMC) can be deformed and deep into the lower half complex plane with the endpoints q = 0, Λ fixed.And then the deformed path of E ′ can avoid crossing the branch cut parallel to the negative imaginary axis, given by F (see Fig. 7 in Ref. [100]).Another advantage of such a parametrization is that the divergence of the pion-exchange potentials can also be moderated to some extent [128] which we will discuss in the next subsection.

The particle-exchange potentials
Due to the unitarity, the t and u-channel particle-exchange potentials can always lead to the LHCs when projecting these potentials into partial waves.For a general process a(p 1 )b(p 2 ) → c(p 3 )d(p 4 ), the t-channel S-wave projected potential is, with m i /E i /p i (i = 1, 2, 3, 4) the mass/energy/module of the three-vector momentum of particle a, b, c, d in order, and m ex the mass of the exchanged particle.Note that the last second step to the last one is reasonable only for p i on the physical axis or around it.Otherwise, the potential in such a form has to be analytically extrapolated properly.The branch cuts are thus determined by restricting the arguments of the two logarithms into (−∞, 0), i.e., where the two "±" signs are uncorrelated.Equation (B5) gives rise to the curved cut, called dynamical cuts (DCs), with p 3 as a function of p 1 and the exact discontinuities can be evaluated by the N/D method [127].As for t-channel , the on-shell, half-off-shell and off-shell projected partial potentials are argued as follows: 1.For the on-shell case p 1 ≈ p 3 , Eq. (B5) with minus sign in the first ± symbol gives, 2. For the half-off-shell case p 1 > 0, the imaginary part of p 3 , dominated by m ex , is still far from the physical region of interests (about tens of MeV mostly) since the masses of the t-channel exchanged mesons are usually at the amount of hundreds of MeV (ρ, ω, σ, etc).
3. For the off-shell case with p 1 , p 3 > 0, there is no crossing of cuts as seen from the second line in Eq. (B4).
While for the u-channel processes (by performing p 3 ↔ p 4 ), the situation turns to be subtle as |E 4 − E 1 | becomes comparable with or even larger than m ex for the exchanged pion while still the same for heavier exchanged mesons.Similarly, the projected partial potentials are treated as follows: 1.For the on-shell case, the poles of the integrand ( 2. For the off-shell case with p 1 , p 4 > 0, the above condition turns out to be a finite cut (0, ] in the √ s complex plane (which fully covers (0, m D 0 + m D * + ] for process D 0 D * + π + − − → D * + D 0 ) and thus crossing over the cut will occur which produces discontinuity when the energy √ s moves from the upper half plane to the lower one.Thus, one needs to deform the integral path rather than [−1, 1] on the real axis such as a polyline −1 → −1 − ia → 1 − ia → 1 with sufficiently large a as used in Ref. [95] or analytically continue by adding the discontinuity on the real axis directly, i.e., twice of the imaginary part on the cut, 3. The situation becomes even more complicated for the half-off-shell case.But the procedure goes the same as the previous case with the mentioned adjustment above.
In our work, we have made a proper treatment for the analytic issues as discussed above, and in Fig. 10 the singular region of the integrand for l ∈ Γ SMC and z ∈ [−1, 1] is illustrated.We have to admit that there are still certain flaws in our numerical treatment though they do not affect what we have achieved.It is worth noting that a better treatment of the on-shell OPE is to decompose it into two parts in TOPT as in Fig. 2(b) since the OPE receives its cut only from the forward diagrams, i.e. the DDπ cuts in our work [96].In this subsection, we give a detailed explanation of the dependence behaviour between Λ and C V found in Fig. 6.To proceed, we simplify the issue of double-charm tetra-quark T + cc into a single-channel scattering problem where the NLO three-body cuts are also turned off.In such a case, the pole of the amplitude, solved by a typical LSE, is evaluated by the root of the determination, namely, 0 = 1 − ∞ 0 d 3 q (2π) 3 V (q)G(q; E)Θ(Λ − q) (C1) = 1 + Λ 0 dq µ DD * q 2 π 2 (2µ DD * E B + q 2 ) V (q) (C2) with µ DD * the reduced mass of the composite DD * system and E B ≈ 273 keV the binding energy of T + cc from the LHCb's Briet-Wigner fitting as an illustration.The typical momentum scale γ T = √ 2µ DD * E B ≈ 22 MeV is far below a typical hard scale 0.3 ∼ Λ ∼ 1.2 GeV.The potentials involved in our work can be parametrized as, and then each integral term in Eq. (C2) is listed as follows, Since the vector-meson-exchange plays the dominant role in the heavy meson-(anti)meson systems, we plot the above three functions and their derivatives with µ i/j ≈ m ρ in Fig. 11.It shows that these three types of potentials lead to different dependence behaviors at different momentum regions: 1) that of contact potential is always proportional to Λ due to the separation Λ ≫ γ T ; 2) that of E1-type potentials, which manifest for a small q, increase with a decreasing rate; 3) that of M1-type potentials, which manifest for a large q, increase with an increasing rate.And then, the dependence in Fig. 6 can be understood similarly if flipping the Λ − C V axes where a stationary point shows up from the competition between f E1 and f M 1 around Λ ≈ 1.0 GeV (the precise position is also related to the prefactors a i /b i /c i ).

FIG. 2 :
FIG. 2: Three-body cuts involved in the DD * scattering.(a) Self energy of D * by the Dπ loop; (b) On-shell decay of process D * → Dπ.The left hand side is the propagator of pion in Feynmann representation used in our work and the right hand side is the sum of the forward and backward emissions of pion in TOPT.
is the magnitude of three-momentum in the c.m. frame of D * with E ′ the energy in the c.m. frame of D * and m D and m π the masses of D and pion, respectively, as depicted in Fig.2(a).The first/second Riemann sheets are also defined regarding the imaginary part of p cm induced by the right-hand cut (RHC) lying along [m D + m π ,∞), i.e., RS-I : ℑp cm > 0 , RS-II : ℑp cm < 0.

FIG. 7 :
FIG. 7: Dependence of the binding energy of D * D * (I = 0) on C V (Λ) in the three schemes.The long-distance OPE is removed in the Scheme-I and Scheme-II while retained in the Scheme-III correspondingly.

FIG. 8 :
FIG. 8: Real and imaginary parts of pole positions relative to the D 0 D * 0 + c.c threshold of isoscalar D D * with the positive (upper panels) and negative (lower panels) C-parity.

FIG. 9 :
FIG.9:The real part of p cm of D * + [D 0 π + ]D 0 before(a)/after(b) the analytic continuation within the energy range −500 keV ≤ ℜE ≤ 500 keV relative to D 0 π + D 0 threshold and width range −500 keV ≤ ℑE ≤ 500 keV when the frame D * + [D 0 π + ] is static, i.e., q = 0.The three-body cut is redirected from the RHC on the real axis to the direction parallel to the negative imaginary axis.

) which leads to the LHC along p 2 3 ≤ − m 2 ex 4 ,
which means that the crossing of cuts only happens when |ℑp 3 | ≥ mex 2 .
0 for z within the interval [−1, 1] actually lead to different branch cuts for the three different processes involved: 1) a finite cut on the real axis above the D 0 D * + threshold for process D 0 D * + π + − − → D * + D 0 ; 2) a finite cut along with a circular cut as shown in Fig. 22 in Ref. [100] for process D 0 D * + π 0 −→ D * 0 D + ; 3) a finite cut crossing the D 0 D * + threshold for D + D * 0 π + − − → D * 0 D + which can be treated analytically as shown later.

FIG. 10 :
FIG.10:The complex contour Γ SMC (red line) and its singular region (gaussian points in light blue or purple) given by Eq. (B8) for processD 0 D * + π 0 −→ D * 0 D + (left) and D 0 D * + π − − − → D * + D 0 (right).The parameters used in this demo are Λ = 0.5 GeV, ω = 0.1 GeV and V 0 = −0.1 GeV.As shown in the figures, there is no analytical problem with the left process when solving the LSE.However, crossing the cut [−1, 1] will occur for the right one, which can be treated properly by Eq. (B7).See the text for further explanations.
Appendix C: Correlation between cutoff parameter Λ and potential strength CV

15 FIG. 11 :
FIG.11:The functions (left) and their derivatives (right) of the three types of potentials convoluted with a two-body Green function.

TABLE I :
Parameters and pole positions of the T + cc relative to D 0 D * + threshold.The first Riemann sheet (RS