Heavy quark spin partners of the $Y(4260)$ in coupled-channel formalism

The charmoniumlike state $Y(4260)$ is described as predominantly a $D_1 \bar{D}$ molecule in a coupled-channel quark model [Phys.\,Rev.\,D\,\textbf{96},\,114022\,(2017)]. The heavy quark spin symmetry (HQSS) thus implies the possible emergence of its heavy quark spin partners with molecular configuration as $D_1 \bar{D}^*$ and $D_2^* \bar{D}^*$ below these charmed mesons' thresholds. We analyze the probabilities of various intermediate charmed meson loops for $J^{PC}=1^{--}$ exotic state $Y(4360)$ and find that the channel $D_1 \bar{D}^*$ couples more strongly around its mass regime, and the coupling behavior remains the same even if the mass of $Y(4360)$ is pushed closer to $D_1 \bar{D}^*$ threshold. This enlightens that the most favorable molecular scenario for the $Y(4360)$ could be $D_1 \bar{D}^*$, and hence it can be interpreted as HQSS partner of the $Y(4260)$. We also find the strong coupling behavior of $D_2^* \bar{D}^*$ channel with the $\psi(4415)$, which makes it a good candidate for a dominant $D_2^* \bar{D}^*$ molecule. We discuss the important decay patterns of these resonances to disentangle their long- and short-distance structures.


I. INTRODUCTION
Heavy quark spin symmetry (HQSS) implies that the hadronic interactions do not depend on the spin of the heavy quark, and heavy mesons can be fully classified by using the quantum numbers of the light quark cloud. If a hadron with a given heavy quark spin is observed experimentally, it is then natural to expect that there should exist its spin partners with different heavy quark spin but with the same light degree of freedom. Such heavy quark spin patterns are identified for instance for recently observed LHCb pentaquarks [1][2][3][4][5][6], and for charged exotic states in the bottom sector, namely Z b s [7][8][9][10][11][12].
In this work, we predict the heavy quark spin partners of the exotic state Y(4260), provided that this state is predominantly a D 1D hadronic molecule 1 [13]. The HQSS thus implies the possible emergence of other hadronic molecules such as D 1D * and D * 2D * below or nearby these thresholds. We demonstrate that, as Y(4260) is dominant by the D 1D molecular component, its HQSS partners D 1D * and D * 2D * then naturally emerge and we identify them as J PC = 1 −− charmoniumlike state Y(4360) and ψ(4415), respectively. It is wort noting that the mass difference between Y(4360) and Y(4260) is almost same as the mass difference between vector and pseudoscalar charmed meson, namely This makes a good benchmark to expect Y(4360) as an ideal candidate for the HQSS partner of the Y(4260). The nearest threshold involving two charmed mesons above Y(4360) is D 1D * (see Fig. 1 for other threshold levels), and it is approximately as far as D 1D from the Y(4260) 2 . This reflects an important consequence of the HQSS-the binding energies of these meson molecules would be of the same order [15]. Therefore, one can expect that the coupling strength of Y(4360) with the first few channels will be of the same order as we observed for the Y(4260) case. We will briefly come back to this point in Sec. III. Moreover, HQSS implies the presence of D 1D and D * 2D components in the wave function of Y(4360), the latter channel couples to J PC = 1 −− through D-wave only. However, their masses are below the experimental value of the mass of Y(4360), and hence will be neglected here. The mass values of relevant thresholds can be read from the Appendix (Table VI).
Before we go into details of our calculations, we first briefly review the status of Y(4360). This charmoniumlike structure was first observed in 2007 at BaBar in the initial state radiation (ISR) process e + e − → π + π − ψ(2S ) [16], and later, in the same year at Belle [17]. Subsequently, BaBar and Belle have updated their data [18,19]. Recently, BESIII Collaboration has observed this state for the first time in the e + e − → π + π − J/ψ process [20]. We summarize the extractions of mass and width of the Y(4360) from different experimental data in Table I. More details on its status can be found in a recent review [21], and for different theoretical interpretations, see Sec. 4.8 of the review article [22]. Y(4360) is recognized as a potential candidate for an exotic state since it does not show strong coupling to open-charm channels (such as DD) which are generally expected as dominant decays of vector charmonia. Moreover, there is not any pronounced enhancement around the Y(4360) mass in the in-    [20] is based on the e + e − → π + π − J/ψ process.
clusive cross sections e + e − → hadrons, the so-called R-value measurements. Hence, it is necessary to make more effort to investigate the structure of Y(4360) with updated knowledge of multiquark dynamics. We use this opportunity and report our analysis of this exotic state. The J PC = 1 −− exotic state Y(4260) is being analyzed in the coupled-channel quark model [24], where the constituent quark model is used to describe the bare quark-antiquark interaction, and the 3 P 0 quark-pair creation mechanism is used to couple charmonium core to the molecular components. The probabilities of the various charmed meson components of the Y(4260) are analyzed and it is found that, even though the HQSS forbids S -wave coupling of D 1D to the 3 S 1 charmonia [ψ(nS )], the D-wave coupling is allowed and not negligible. The Y(4260) is interpreted as a mixture of a charmonium core plus the dominant D 1D component. To probe the charmonium core of Y(4260), the coupling behavior of the D 1D channel with different charmonium cores is investigated and it is found that it couples more strongly to 3 D 1 charmonia [ψ(nD)]. This manifests that the charmonium core of the Y(4260) is likely to be ψ(nD). The production of ψ(nD) via e + e − annihilation is suppressed, that explains why the production cross section of e + e − → vector charmonia exhibits a dip around 4.26 GeV.
In this work, we extend our coupled-channel formalism to investigate other J PC = 1 −− exotic states. The purpose is to build up a heavy quark spin multiplet of charmed meson molecules [25]. In what follows, we do not aim to fully explain the production or decay patterns of these exotic hadrons. Instead, we will investigate the molecular components in experimentally observed J PC = 1 −− exotic states. In unquenched quark model (UQM), the wave functions encapsulate both long-and short-distance information, and we intend to utilize them to explore the structure of exotic Y and ψ resonances.
A brief introduction of our coupled-channel framework is given in Sec. II, where some subtleties of the HQSS are also discussed. Section III is devoted to the analysis of our results and their interpretations. In Sec. IV, a few remarks on the strong decays of Y(4360) and ψ(4415) are provided, followed by a summary of this study in the Sec. V.

II. COUPLED-CHANNEL FORMALISM
The quenched quark model described the spectrum of the low-lying states reasonably well, but its predictions for the states nearby or above open-flavor thresholds are questionable. It lacks the influence of "sea-quarks". Heavy quarkonium can couple to intermediate heavy mesons through the creation of light quark-antiquark pair. This enlarges the Fock space of heavy quarkonium and manifests the presence of multiquark components in its wave function. Such components will change the Hamiltonian of the potential model, causing a mass shift due to self-energy corrections, and may also give direct contributions to strong and electromagnetic decays [26,27]. The probability of such multiquark (molecular) components can be worked out (see e.g., Ref. [24]).
In UQM, a physical or experimentally observed hadron |A can be expressed as where c 0 and c BC stand for the normalization constant of the bare state and the BC components, respectively. In this work, B and C refer to charmed and anti-charmed mesons, and the summation over BC is carried out up to ground state P-wave charmed mesons [28]. The bare state |ψ 0 is normalized to 1, and the physical state |A is also normalized to 1 if it lies below DD threshold. |BC; p is normalized as BC; The effects from the BC components are referred to as coupled-channel effects. The full Hamiltonian of the physical state reads as where H 0 is the Hamiltonian of the bare state (see Appendix A for details), the continuum Hamiltonian is H BC |BC; p = E BC |BC; p with E BC = m 2 B + p 2 + m 2 C + p 2 is the energy of the continuum state (the interaction between B and C mesons is neglected here and transition between one continuum to another is not included), and H I is the interaction Hamiltonian which triggers the mixing of the bare state to the continuum.
For the bare-continuum mixing, which is an important dynamical pieces of the UQM, we adopt the widely used 3 P 0 model [29]. In this model, the nonperturbative creation of light quark-antiquark pairs is triggered from the vacuum of quantum chromodynamics (QCD) with J PC = 0 ++ , which in the spectroscopical notation 2S +1 L J can be written as 3 P 0 [26,30]. The interaction Hamiltonian can be expressed as where m q is the produced quark mass, and γ is the dimensionless coupling constant. The ψ q (ψ q ) is the spinor field to generate anti-quark (quark). Since the probability to generate heavier quarks is suppressed, we use the effective strength γ s = m q m s γ in the following calculation, where m q = m u = m d is the constituent quark mass of up (or down) quark and m s is strange quark mass. Their numerical values are listed in Appendix A (Table V).
The mass shift caused by the BC components and their probabilities are obtained after solving the Schrödinger equation with the full Hamiltonian H. They are expressed as where M and M 0 are the eigenvalues of the full (H) and quenched/bare Hamiltonian (H 0 ), respectively; P BC is the unnormalized probability, which is also called the coupling strength in the next section. In order to analyze different partial-wave contributions, we adopt the Jacob-Wick formula to separate different partial waves of P BC [31]. (See Ref. [26] for a derivation of the above relations and UQM calculation details.) To proceed, we need to specify that the coupled-channel calculation cannot be pursued if the wave functions of the |ψ 0 and BC components are not settled in Eqs. (5)−(6). The major part of the coupled-channel calculation is encoded in the wave function overlap integration, , and m Q and m q denote the charm quark and the light quark mass, respectively. The φ 0 , φ B and φ C are the wave functions of |ψ 0 and BC components, respectively and the notation * stands for the complex conjugate. These wave functions are in momentum space, and they are obtained by the Fourier transformation of the eigenfunctions of the bare Hamiltonian H 0 . For the heavy quarkonium decays, the heavy quarks are treated as spectators in the 3 P 0 model. Their polarizations will not change after the generation of the light quark-antiquark pairs. This leads us to conclude that the 3 P 0 model itself respects the HQSS. Nevertheless, some HQSS breaking effects can still slip into the coupled-channel calculation. The breaking effects mainly lie in the input of the charmed mesons masses and can be noticed from their mass splitting in the same j P l multiplet, where j l is the total spin of the light-quark cloud (with parity P) in the charmed meson. For example, D * 2 and D 1 belong to the same j + l = 3/2 multiplet, but they do not degenerate experimentally as claimed by the HQSS.
The HQSS further leads to configuration mixings of charmed mesons having different j l . For example, the experimentally observed D 1 (2420) (or D 1 in short) and D 1 (2430) (or D 1 ) are not the pure j + l = 3/2 and j + l = 1/2 (or the quark model 3 P 1 and 1 P 1 ) states, respectively, but are their linear combinations. Therefore, the wave functions from the quark model should be modified accordingly. The mixture can be formulated as where θ is the mixing angle. The HQSS predicts it to be θ 0 = arctan( √ 2) ≈ 54.7 • , which is called the ideal mixing angle [32,33]. Experimentally, the deviation from ideal mixing is observed (5.7 • ± 4 • ) [34], which is very small compared with θ 0 ≈ 54.7 • . As reported earlier [24], our results have minor changes if this deviation is added to θ 0 . Hence, for simplicity, we will use ideal mixing of HQSS in this study.
Our bare mass value for the ψ(3D) is nearly 250 MeV higher than the Y(4360), as listed in Table V. Hadron loops can shift bare mass down to Y(4360) to compensate such a large mass gap [35]. This requires the fitting of the charmonium spectrum which involves much ab initio calculations (such as fixing of 3 P 0 's coupling constant γ) which is beyond the scope of this work. Therefore, γ is not fixed here, that explains why P BC in Eq. (6) is unnormalized. As a consequence, the absolute probabilities of the various charmed meson components are not determined. Instead, by analyzing the P BC , one can still draw very useful conclusions [24], as demonstrated in the following sections.  [36,37]. To compare the results under different assumptions, we calculate the probabilities for all four aforementioned bare states. For consistency, the same set of parameters is being adopted which was previously used in the study of Y(4260) [24].
The unnormalized coupling strengths of Y(4360) with the charmed meson components are given in Table II. Due to higher mass of Y(4360), we only consider those channels which involve at least one P-wave charmed meson. It is interesting to notice that one of the channels with largest coupling to Y(4360) is D 1D * . The channel D * 2D * also shows a sizable coupling with Y(4360) for ψ(3S ) and ψ(4S ) cores. However, due to several reasons, discussed in the following, Y(4360) can be regarded as having a dominant D 1D * component. For the coupled channels which involve one S -and one P-wave charmed meson, the HQSS leads to set . As a result, the contribution of all those channels which involve charmed meson form the same j l multiplet will be of the same order [except for the D − D 1 (D 1 ) and D * 2D channels, which are the decay channels of Y(4360)]. In such a scenario, the largest coupling will naturally come from the D * 2D * channel, since the spin configurations for this channel are at the maximum. We refer this as spin enhancement mechanism. However, if the physical masses are applied to these charmed mesons, the D * 2D * channel will be a little further from the Y(4360) and is expected to have a smaller coupling strength.
For Y(4360), we noticed that the spin-enhancement mechanism is giving sizable contribution when we use the experimental mass of Y(4360) as the PDG average value 4368 ± 13 MeV. However, a recent analysis by the BESIII Collaboration for the process e + e − → π + π − ψ(2S ) reported a larger mass of this resonance as 4383.8 ± 4.2 ± 0.8 [23], and an independent analysis of the BESIII combined data of e + e − → π + π − J/ψ and e + e − → π + π − ψ(2S ) extracted the mass of Y(4360) resonance as 4386.4 ± 2.1 ± 6.4 [38]. This indicates that Y(4360) might have a larger mass. The coupling of Y(4360) with different charmed meson channels found to be sensitive to its mass. With the use of a larger mass, the coupling of Y(4360) to the channel D 1D * dominantly exceeds the D * 2D * channel and makes it a more reasonable candidate for a D 1D * molecule. The couplings to those channels that are far above 4.368 GeV are generally small-a universal conclusion from the coupled-channel effects. The asymptotic behavior of P BC is proportional to 1/(m B + m C ) 2 . If the coupled channels are further from Y(4360), their contributions will be naturally suppressed; we call this the mass-suppression mechanism. The emergence of this mechanism can be seen from Table II, where the higher channels have very small contributions.
For a solid conclusion, the two mechanisms-mass suppression and the spin enhancement-have to compete with each other to tell us which channel gives the dominant contribution. We will come to this point in the next subsection. If the D s mesons are involved in the coupling channels, then an additional suppression comes from the effective strength of the 3 P 0 model γ s . Since γ s ≈ 0.66γ (using the constituent quark mass from Appendix A), the couplings to D s mesons are universally smaller than the D mesons.
It is worth mentioning that, even though the non-D 1D * components are suppressed when the charmonium core is ψ(nD), the contribution of these components could still be sizable. As shown in Table II, for the ψ(3D), the contribution from the D * 2D * channel is still around the D 1D * one. However, the mass prediction for the quark model ψ(3D) state is nearly 200 MeV − 300 MeV above the Y(4360) [36]; thus, having the  ψ(3D) core in its wave function is not well justified. The most likely possibility for the charmonium core of the Y(4360) is ψ(2D) due to the following reasons: 1. The production cross section of vector charmonia in e + e − → hadrons (the R-value) is proportional the wave function at the origin, which is zero for the case of ψ(nD). Since the ψ(nD) can only couple to the virtual photon at the next-to-next-to leading order [39], its direct production at the e + e − collider is suppressed.
2. The R-value measured by the BES Collaboration [40] has a dip instead of a peak around 4.368 GeV, which indicates that the Y(4360) is likely to have a D-wave charmonium core, and its mass value makes ψ(2D) the most promising.
3. The coupling of ψ(2D) with the D 1D * is found to be the maximum, in two different sets of parameters and wave functions (see details in Table III). In this sense, our results support the dominant long-distance component of Y(4360) to be the D 1D * .
This enables us to conclude that the structure of Y(4360) is a mixture of a charmonium core ψ(2D) with the dominant molecular component D 1D * . To cross check the model dependence of our results and the validity of our conclusions, we try different sets of parameters and different wave-function approximations. For this purpose, we adapt the quark model parameters from Barnes and Swanson [35] and simple harmonic oscillator (SHO) approximations for the wave functions. The results for the normalized coupling strengths are compared in Table III. The wave functions used by Barnes and Swanson [35] are SHO approximations [41] but are useful to cross check the coupling pattern. One of the largest coupling channels to Y(4360) is D 1D * even if we use different set of parameters. This conclusion, to some extent, is model independent, since we tried two different sets of parameters and different wave functions to cross check the model dependence.
The behavior of the D 1D * coupling strength as a function of mass of the Y(4360) resonance is shown in Fig. 2. We try different initial charmonium wave functions and push the initial mass of the resonance to vary in a vast energy range up to 4.4 GeV. A stable behavior of the coupling is found. The coupling gets enhanced as the initial mass of the resonance approaches the D 1D * threshold, and this behavior is the same for all considered initial charmonium wave functions. This verifies our conjecture, i.e., Y(4360) has the largest coupling to the D 1D * channel in all wave-function choices. Hence, we can conclude that Y(4360) can be described in a molecular picture having D 1D * as its largest component, and thus can be interpreted as a spin partner of the Y(4260).

2D * Channel
For the case of ψ(4415), since it is extremely close to the D 1D * threshold, one can expect that this state will show dominant coupling to the aforementioned channel. However, it is found that the spin-enhancement mechanism takes over here and the ψ(4415) has shown the largest coupling to the D * 2D * channel (see Table IV). The coupling of this channel with the ψ(4S ) and ψ(2D) cores is almost the same. The behavior of the D * 2D * coupling strength as a function of mass of the ψ(4415) resonance is shown in Fig. 3.
Since the next higher channels above D * 2D * , such as D * 2D * 0 or D * 2D * 2 (both of these channels couple to J PC = 1 −− in Pwave), are further from the ψ(4415), and the mass suppression mechanism implies that the contributions from all these higher components will be highly suppressed, which makes the name of the D * quark spin structure such as D 1D and D * 2D ; the latter channel couples to J PC = 1 −− through D-wave only. However, these channels are open for ψ(4415) because the mass of these thresholds are below the experimental mass of ψ(4415). This is merely the reason for zero coupling strengths in the Table IV.
An important observation is that, if the charmonium core of the ψ(4415) is in the S -wave, the two mechanisms (mass suppression and the spin enhancement) will be of similar importance. As a consequence, non-negligible contributions other than the D * 2D * channel may exist. This indeed showed up in our results and can be seen from Table IV. The coupling strength of ψ(4415) to the D 1D * is very close to the D * 2D * channel for the case of ψ(3S ) and ψ(4S ). In such a scenario, one can argue that the ψ(4415) is dominant by D 1D * rather than D * 2D * . However, the sizable coupling of the ψ(4415) with the channel D 1D * is merely an artifact of just being extremely close to this threshold.
It is most likely that the charmonium core of ψ(4415) is in S -wave which results the same order of coupling for all initial core wave functions, as can be seen from Table IV. This is further supported by the R-value measurement that the production cross section is enhanced significantly around ψ(4415)'s mass which normally is the case for ψ(nS ) cores. The mass difference between ψ(4415) and the D * 2D * threshold is just close enough to have the same binding energy as Y(4260) which reflects an important consequence of HQSS.
This analysis leads us to conclude that ψ(4415) is likely to be a mixture of a short-distance core ψ(nS ) and the dominant long-distance D * 2D * component. However, coupling strengths alone are not enough to conclude the structure of these resonances. An important way to disentangle the molecular configuration of Y(4360) and ψ(4415) is to look for suggested decay patterns, as discussed in the following section.

Coupled Channels
Benchmark-I Benchmark-II   Channels

IV. REMARKS ON DECAYS
The open-charm decays of Y(4360) provide an important pathway to disentangle its long-distance structure. The natural decay of a molecular state is into its constituents [15], and the strongest decay channel of the D 1 is D * π. Hence, in the D 1D * molecular configuration, we argue that the dominant decay mode of the Y(4360) is likely to be the D * D * π, and the decay of Y(4360) into the DD * π final state is not possible. This is the reason that the Y(4360) has not shown up in the e + e − → DD * π final state at Belle [42].
In fact, the DD * π final state is expected to be the dominant decay mode of the Y(4260) in the molecular picture [43], which is supported by the recent precise measurements of BESIII Collaboration [44]. It is interesting to mention that the Y(4360) is again invisible in this recent BESIII data [44]. Hence, it is very important to measure e + e − → D * D * π cross sections to prove/exclude the existence (as a hadronic molecule) of Y(4360) which we claim is a heavy quark spin partner of the Y(4260).
In contrast to above, the D * 2 meson can decay into D * π and Dπ. Therefore, a D * 2D * molecular state must leave strong imprints in the DD * π and D * D * π final states. For ψ(4415) → DD * π, an upper limit for the partial decay width has been extracted by the Belle to be smaller than 11% [14,42]. The broad enhancement around 4.40 GeV in the recent BESIII measurement of e + e − → DD * π [44] indicates clear contributions from ψ(4415). If the same pattern is observed in the e + e − → D * D * π process, it will indicate that the ψ(4415) is likely to have a dominant D * 2D * component which would provide a quantitative support to our conclusion.
The precise evaluation of the decay widths of the above discussed decays requires an accurate knowledge of the coupling of a molecular state with its components. One can benefit from the recent BESIII data [44] to extract the Y(4260)−D 1D effective coupling and use HQSS to relate it to other molecular configurations. This would be the topic of our future explorations [45].

V. SUMMARY
We computed and analyzed the probabilities of various charmed meson molecular components for Y(4360) under the coupled-channel formalism by assigning ψ(3S ), ψ(4S ), ψ(2D) and ψ(3D) initial wave functions. We found that the channel D 1D * couples more strongly around the Y(4360) mass regime, and the coupling behavior is the same for all considered initial charmonium wave functions. This enlightens that the most favorable molecular scenario for Y(4360) is the D 1D * , and hence, it can be interpreted as a heavy quark spin partner of the Y(4260).
We also analyze the coupling of ψ(4415) with several nearest charmed meson channels. It shows the largest coupling to the D * 2D * channel in all four initial charmonium wave functions, which makes ψ(4415) a good candidate for a prominent D * 2D * molecule. By this means, we argue that the heavy quark spin multiplet involving one P-and one S -wave charmed meson is considered as completed.
However, the hitherto unobserved decays of these resonances are highly decisive and demanding, such as Ψ → D * D * π, where Ψ ∈ {Y(4360), ψ(4415)}. Once the predicted pattern of heavy quark spin multiplet is confirmed by future experiments, it will serve to deepen our understanding of how QCD forms hadronic matter by arranging multiquarks.   [14].