Light- and heavy-quark symmetries and the $Y(4230)$, $Y(4360)$, $Y(4500)$, $Y(4620)$ and $X(4630)$ resonances

The heavy hadron spectrum is constrained by symmetries, of which two of the most important ones are heavy-quark spin and SU(3)-flavor symmetries. Here we argue that in the molecular picture the $Y(4230)$ (or $Y(4260)$), the $Y(4360)$ and the recently discovered $Y(4500)$ and $Y(4620)$ vector-like resonances are linked by these two symmetries. By formulating a contact-range effective field theory for the $D \bar{D}_1$ and $D_s \bar{D}_{s1}$ family of S- and P-wave charmed meson-antimeson systems, we find that if the $Y(4230)$ were to be a pure $D \bar{D}_1$ molecular state, there would be a $D^* \bar{D}_1$ partner with a mass similar to the $Y(4360)$, a $D_s \bar{D}_{s1}$ partner with a mass close to the $Y(4500)$ and three $J=1,2$ $D_s^* \bar{D}_{s1}$ and $J=3$ $D_s^* \bar{D}_{s2}^{*}$ bound states with a mass in the vicinity of $4630\,{\rm MeV}$, of which the first one ($J=1$) might correspond with the $Y(4620)$. The previous predictions can in turn be improved by modifying the assumptions we have used to build the effective field theory. In particular, if we consider the closeness of the $D^* \bar{D}_1$-$D^* \bar{D}_2^*$ and $D_s^* \bar{D}_{s1}$-$D_s^* \bar{D}_{s2}^*$ thresholds and include the related coupled channel dynamics, we predict a $J=2$ positive C-parity state with a mass around $4650\,{\rm MeV}$. This hidden-strange and hidden-charm state might in turn be identified with the $X(4630)$ that has been discovered past year by the LHCb in the $J/\psi \phi$ invariant mass distribution.


I. INTRODUCTION
The Y (4260) -or Y (4230), as later measurements have converged to a lower value of its mass -and Y (4360) were originally observed by the BaBar collaboration in the initial state radiation process e + e − → γ ISR J/ψπ + π − [1] and in the e + e − → ψ(2S)π + π − invariant mass distribution [2], respectively. Both are charmonium-like vector states that were not expected in the quark-model and their nature still remains elusive. The Y (4230) has been proposed to have a sizable D 1D molecular component [3][4][5][6][7][8][9], which is the premise we will use in this work, though there are other interesting conjectures about its nature [10][11][12][13][14][15][16][17]. The Y (4360) is also suspected to have non-charmonium components [13, 15,18]. A recent observation in ηJ/ψ by the BESIII collaboration [19] provides a mass and width of for the Y (4230) and for the Y (4360), where we notice that in [19] these two resonances are referred to as Y (4220) and Y (4390), following the nomenclature of their previous observation in π + π − h c [20]. However, the experimental situation is not * yanmaojun@itp.ac.cn † mpavon@buaa.edu.cn necessarily clear with respect to whether the Y (4230) and Y (4360) are each a unique states or not, with the 2020 edition of the Review of Particle Physics (RPP) [21] subdividing each one into two entries, ψ(4230), ψ(4260) for the Y (4230) and ψ(4360), ψ(4390) for the Y (4360) (while the more recent 2022 edition does not [22]). Here we will treat the Y (4230) and Y (4360) as if each were a unique state.
It is interesting to notice the most recent observation of the Y (4230) by the BESIII collaboration [23] (in e + e − → K + K − J/ψ), which gives Γ Y = (72.9 ± 6.1 ± 30.8) MeV , and finds a new resonance with a mass and width of M Ys = (4484.7 ± 13.3 ± 24.1) MeV , which is referred to as the Y (4500), in reference to its mass. In this work we will conjecture that the Y (4500) is the hidden-strange partner of the Y (4230). Yet, there are other two hidden-charm resonances of interest to us. The Belle collaboration has observed [24] a resonance resembling the Y (4230) in the e + e − → D + s D s1 (2526) − + c.c process, with mass and width M Y * s = (4625.9 +6.2 −6.0 ± 0.3) MeV , Γ Y * s = (49.8 +13.9 −11.5 ± 4.0) MeV , which we will refer to as Y (4620). This discovery has been followed by a very similar resonance also found by Belle [25] in the e + e − → D + s D s2 (2573) − + c.c process with compatible mass and width which we will presume to be the same state as the one found in [24]. Finally, the LHCb collaboration has discovered a series of ccss and ccus resonances in the J/ψφ and J/ψK + invariant mass distributions [26], of which we will highlight the X(4630): This discovery connects with previous predictions of D * sD * s2 virtual states with a similar mass [27], but it has also been explained as a bound state [28] and there are non-molecular explanations too [29,30].
Here we will explore the possible relation among these five states within the molecular picture, where we will assume that they are S-and P-wave charmed mesonantimesons bound states (where, for the P-wave charmed mesons, we will be referring to the ones with total lightquark spin J L = 3 2 ). The tentative molecular interpretations we will propose for these five states are listed in Table I, together with the binding energies. For describing them as molecules we will formulate a non-relativistic effective field theory (EFT) which incorporates SU(3)flavor and heavy-quark spin symmetry (HQSS), in line with previous EFT works for S-wave charmed mesonantimeson systems [31][32][33][34][35][36][37].
Before presenting the complete derivation of our formalism, we will provide a brief, heuristic overview of our results. If we consider the Y (4230) and Y (4360) as 1 −− DD 1 and D * D 1 molecular states and in a first approximation ignore the kinetic energy of the charmed mesons (as they are heavy), their masses will be given by m(Y (4360)) = m(D * ) + m(D 1 ) − V 0 (Y * ) , (16) where V I (Y ( * ) ) is the expected value of the potential binding the previous meson-antimeson systems in the isospin I = 0 configuration. Conversely, if the Y (4500) and Y (4620) happen to be 1 −− D sDs1 and D * sD s1 molecules, their potentials should be directly related to those of the Y (4230) and Y (4360) via HQSS and SU(3)-flavor symmetry. If we additionally assume that the potentials in the I = 1 configurations are negligible in comparison with the I = 0 ones (e.g. if the potential is given by vector-meson exchange), we end up with m(Y (4500)) = m(D s ) + m(D s1 ) − 1 2 V 0 (Y ) , (17) m(Y (4620)) = m(D * s ) + m(D s1 ) − This means that we can calculate the masses of the Y (4500) and Y (4620) resonances from the ones of the which are in line with their experimental masses, see Eqs. (7) and (11). This suggests that the Y (4500) and Y (4620) are simply the hidden-strange partners of the Y (4230) and Y (4360). In contrast, the explanation of the X(4630) is more involved and we will discuss it later, as it requires to consider the explicit HQSS structure of the potential as well as the inclusion of coupled channel effects.
In the following lines we will explain in more detail the previous arguments and the constraints imposed by HQSS and SU(3)-flavor symmetry to the family of molecular states composed of an S-wave and P-wave charmed meson-antimeson pair. The manuscript is structured as follows: in Sect. II we write down the lowest order contact-range potential compatible with HQSS for these systems; in Sect. III we consider the effects of SU(3)flavor; in Sect. IV we propose two possible power countings by which to organize our limited knowledge of these systems. Finally, in Sect. V we summarize our findings. In addition Appendix A analyzes whether pion exchanges are perturbative, while in Appendix B we explain a few technical details about the spin operators appearing in these bound states and Appendix C describes two phenomenological models to which we compare our EFT results.

II. HEAVY-QUARK SPIN SYMMETRY AND THE CONTACT-RANGE POTENTIAL
We begin by explaining how HQSS constrains the Swave contact-range interaction between the S-and Pwave charmed mesons. The interaction (derived for the first time in Ref. [38]) contains four components and four unknown coupling constants, which increase to eight once we consider isospin and SU(3)-flavor symmetries.
The choice of a contact-range potential is grounded on the observation that this is expected to be the lowest order description within an effective field theory (EFT) with charmed mesons and pions as degrees of freedom [32][33][34].
Pions are perturbative at low energies, as we explain in detail in Appendix A where we compare the sizes of iterated and tree-level one pion exchange (OPE). There we find that the inclusion of non-perturbative OPE effects is not required for momenta below 800 MeV (at the minimum), which basically is of the same order of the ρ mass, i.e. the momentum scale at which we expect the EFT description not to be valid anymore. It is worth noticing that this is analogous to what happens with OPE in two-body systems composed of S-wave charmed mesons [32,34], and compatible with recent calculations of the effect of pions on the binding energy of S-and P-wave charmed meson-antimeson molecules [9], which happens to be small. Consequently, we will consider pions to be subleading.
In contrast, contact-range couplings are leading: the EFT description of shallow bound states requires the promotion of this type of couplings to leading order [39]. Of course there will be limitations with this description: the S-and P-wave charmed meson-antimeson molecular candidates are not shallow, but show binding energies for which the EFT expansion is not expected to converge particularly well. For instance, if we consider the Y (4230) to be a DD 1 bound state then it will have a binding energy of about B ∼ 70 MeV and a binding momentum γ = √ 2µB ∼ 390 MeV (where µ is the reduced mass of the system). From this we expect the expansion parameter of the EFT to be the ratio of this binding momentum over the rho meson mass or γ/m ρ ∼ 0.5. This ratio technically lies within the range of validity of a possible EFT description, but will result in poor convergence.
Considerations of convergence aside, the four components of the S-wave contact-range potential can be characterized as charge, dipolar magnetic, dipolar electric and quadrupolar magnetic in an analogy with the multipolar expansion of electromagnetic interactions (check, e.g. Refs. [40,41]). We first show how this decomposition arises by writing down the meson-meson interaction in the standard superfield notation. Then, by noticing that the interaction only depends on the light-spin of the heavy mesons, we will rewrite the meson-meson contactrange potential in a more compact form using a notation that exploits this observation.

A. Heavy superfield notation
From HQSS we expect heavy hadron interactions to be independent of the spin of the heavy-quarks within them. This symmetry is implemented at the practical level by writing down superfields that group together the ground and excited states of a heavy hadron. For the S-wave charmed mesons D and D * , their superfield is while for the P-wave charmed mesons D 1 and D * 2 (that is, we are considering the P -wave meson for which the total light-spin of the light-quark is J L = 3 2 ) we have where the indices in D 1,j and D * 2,ij refer to their polarization states in the Cartesian basis, while σ i with i = 1, 2, 3 are the Pauli matrices. The superfields H and T i defined here correspond with the non-relativistic limit of the superfields in [42].
With the previous superfields, the most general S-wave contact-range Lagrangian between the H and T fields with no derivatives is and contains four couplings (which become eight once we include isospin or SU(3)-flavor symmetry). For simplicity we have not indicated either whether we are dealing with heavy mesons or antimesons. By expanding this Lagrangian in terms of the heavy meson fields and reordering the indices properly, we will arrive at the potential in Table II. In the particular case of the H and T superfields, this is painstakingly difficult, which is why we will also derive the potential in a second and more convenient notation.

B. Light subfield notation
Now we rewrite the interaction between an S-wave and P -wave heavy meson within the aforementioned second notation. In this alternative notation, which we might call light-quark notation or light subfield notation, we only explicitly consider the spin of the light quarks, while the heavy quarks act basically as spectators. In particular we can represent the S-and P -wave mesons by the non-relativistic light-quark subfields where we indicate the spin and parity of the light-quark degrees of freedom in parentheses. With these fields the lowest-order contact-range Lagrangian between an Swave and P -wave heavy hadron is where C a , C b , D b and D c are coupling constants, σ Li and S Li (with i = 1, 2, 3) the spin operators of the light quarks inside the S-and P-wave heavy mesons (the Pauli and the spin-3 2 matrices, respectively) and Σ L a spin operator involved in the transition from spin-1 2 to -3 2 and vice versa (the matrix elements of which can be consulted in Appendix B), while Q Lij is a tensor spin operator that is defined as: The most important of these couplings will be D b , which we will assume to play a very important role in the description of the Y (4230). From the previous contact-range Lagrangian we obtain the following contact-range potential where we define the product of the two tensor spin operators as This potential is written in terms of the light-spin operators, which do not directly correspond with the heavymeson spin operators. For applying the potential in the heavy-meson basis, we have to supplement the previous potential with a series of rules for translating the lightspin operator into the heavy-meson spin operators. For the S-wave charmed mesons, D and D * , these rules are where J 1 stands for the spin-1 matrices as applied to the D * meson. For the P-wave charmed mesons, D 1 and D * 2 , the rules are where S 1 and S 2 are the spin-1 and -2 matrices as applied to the D 1 and D * 2 mesons, respectively. Finally we have the transition from the S-to P-wave mesons, which involve the Σ L and Q Lij operators. The Σ L matrices are used to represent the D → D 1 , D * → D 1 and D * → D * 2 transitions, which are generated for instance by electric-type dipole E1 vector-meson exchange operators [41]. The transformation rules are plus their hermitian conjugates when we exchange the Sand P-wave mesons, where ǫ is the polarization vector of the D 1 meson, S the spin-1 operator and Σ the operator for the spin-1 to spin-2 transition (which we write down in Appendix B). The choice of different symbols ( J 1 , S 1 and S) for what is essentially the same spin-1 matrices is made to emphasize that they are operating on different states, as it matters where they are sandwiched in between.
Next, the Q Lij matrices represent the D → D 2 , D * → D 1 and D * → D * 2 transitions that can be generated by a magnetic-type quadrupolar M2 vector-meson exchange operator [41]. The transformation rules for the quadrupolar operator Q Lij are where ǫ tij is the tensor polarization of the D * 2 charmed meson, Q the quadrupolar spin operator for spin-1, which is defined as where we remind that S refers to the spin-1 matrices when sandwiched between a D * and a D 1 charmed meson, and Q 12ij is defined as

C. C-parity convention
For calculating the potential in the different J P C configurations it is useful to specify the C-parity convention used, which in our case is where C is the operator transforming particles into antiparticles. That is, we always use the same sign. With Structure of the contact-range potential between a D/D * S-wave charmed meson and a D1/D * 2 P-wave charmed meson. This potential depends on four coupling constants: Ca, C b , D b and Dc. From the saturation arguments of Refs. [40,41], for isoscalar molecules we expect: Ca < 0, C b < 0, D b > 0 and Dc > 0 (check also the discussion in Sect. IV A). If this argument is correct, there will be a series of configurations that are guaranteed to be attractive: this convention, a two-body meson-antimeson state with mesons A and B where A, B are two different combinations of the S-and P-wave charmed mesons considered here, will have Cparity C = η(−1) S−SA−SB , with S the total spin of the AB system and S A , S B the spins of mesons A and B, and η = ±1 a sign.

III. LIGHT SU(3)-FLAVOR SYMMETRY
The light-flavor structure of the contact-range interaction potentially provides further information about heavy meson-antimeson molecules. The constraints that SU(3)flavor symmetry imposes happen to be really simple: as the charmed mesons (antimesons) belong to the3 (3) representation of SU(3)-flavor, their molecules will be the sum of a singlet and octet representation (3 ⊗3 = 1 ⊕ 8). That is, the contact-range potential can be decomposed as a sum of a singlet and octet component, i.e.
where the specific linear combination depends on the isospin and strange content of the molecule. In particular we have where qq and ss refer to the light-quark content of the charmed meson-antimeson system under consideration, with the meaning of q here restricted to q = u, d.
Here two comments are in order: first, the strange and non-strange I = 0 configurations can mix, which means for instance that the Y (4230) and Y (4360) will have a strange component. This type of coupled-channel effect is expected to be small though and we will ignore it. Second, it is worth noticing that the potential in the ss configuration is the average of the isoscalar and isovector ones This happens to be important because from phenomenological models we expect the potential in the isovector configuration to be considerably weaker than in the isoscalar one. The reason is that vector meson exchange cancels out in the I = 1 channel and thus the interaction in this channel depends either on SU(3)-flavor symmetry breaking effects or exchange of more massive mesons (two possibilities being vector charmonia [27,[43][44][45] or axial meson exchange [46]). Naively, we expect these contributions to be weaker than vector meson exchange.

IV. PREDICTIONS AND POWER COUNTINGS
When making predictions with the contact-range theory we have formulated, there will be a significant problem: the lowest-order contact-range potential contains a total of eight unknown couplings (or four if we concentrate in one spin/flavor sector only). In principle the determination of these couplings requires the identification of eight molecular candidates containing different linear combinations of the couplings. Owing to a lack of a suitable number of molecular candidates, this is not a feasible prospect at the moment.
There is a way out of this predicament though: we do not expect all the couplings to be equally important and a few couplings might provide a larger contribution to binding than others. This idea might have support in phenomenology: for instance, vector meson exchange together with vector meson dominance [41] suggest a really strong electric dipole term (and to a lesser extent, magnetic quadrupole term) mediating the D ( * ) → D ( * ) 1 (2) family of transitions. This type of vector meson exchange generates a potential with the same operator structure as the D b term (or the D c term for the magnetic quadrupole transition) in the S-wave contact-range potential. It happens that the D b coupling might be strong enough as to bind the DD 1 system by itself, and the same might be true for the D c coupling and the D * D 1 system. If this is the case, it might be possible to ignore the other couplings at lowest order.
Within the EFT language we do not necessarily worry about the microscopic origin of the contact-range potential we use, be it either light-meson exchange or shorterrange cc components of the wave function [47]. Instead we simply point out that there are several possible power countings, i.e. different ways in which to order the operators of the EFT from more to less relevant. The power counting can be simply assumed, where its suitability and correctness can be later judged on the basis of the predictions a particular power counting entails. In this regard we will explore power countings in which the D b and D c couplings are conjectured to be more relevant than the C a and C b couplings.
A. What we know and what we do not know: proposing a power counting In principle we have eight independent coupling constants, namely C Ia , C Ib , D Ib and D Ic , with isospin I = 0, 1. However, as previously mentioned, there are not enough molecular candidates involving Sand P-wave charmed mesons. The most solid candidate is the Y (4230), theorized to be a 1 −− DD 1 bound state [3][4][5][6][7][8][9] (though its nature is far from clear and other explanations exist [10][11][12][13][14][15][16][17]). With this we can single out one coupling or combination of couplings that is responsible for its binding, and set the other couplings to zero. In terms of power counting this is justified by identifying the chosen coupling as leading order, while all the other couplings will be subleading.
For choosing the coupling, we will first note that the DD 1 potential in the J P C = 1 −− channel contains two couplings where this observation can be combined with phenomenological information. In particular, the S-wave light-meson exchange potential between the S-and Pwave charmed mesons can be written as which at high momenta will receive contributions from the sigma, rho and omega mesons. The explicit form of these components was calculated in [41], yielding (55) where they are written in p-space, with q the exchanged momentum. In the expressions above g Si , g V i and f V i (with i = 1, 2) are the scalar meson and vector meson charge-like and magnetic-like couplings for each of the heavy mesons, f ′ V and h ′ V the electric dipole and magnetic quadrupole couplings of the vector mesons, M a scaling mass (for which a typical choice is the nucleon mass, M = m N = 940 MeV), m S and m V the masses of the scalar and vector mesons and ) the mass difference of the S-and Pwave charmed mesons for the non-diagonal (i.e. W b and W c ) components of the potential. Owing to µ V < m V , the range of the non-diagonal W b and W c components will be larger than that of the the diagonal ones, V a and V b . Also, the V b component differs by a factor of 2/3 to that of Ref. [41] as a consequence of the different normalization used for the spin operators.
At this point it is interesting to notice that the combination of the isospin factor τ 1 · τ 2 and the sign of the omega contribution (ζ = +1 for a heavy meson-meson pair and ζ = −1 for the meson-antimeson case) implies that the vector-meson exchange contributions cancel for I = 1, leading to the conjecture that the potential in the hidden-strange sector is about half that of the isoscalar sector (if we assume that the scalar meson contribution is small), e.g. Eqs. (17) and (18).
Here, before discussing the size of each component, it is important to comment on the dots in Eqs. (52)(53)(54)(55), which we use to indicate Dirac-delta terms that appear as a consequence of the derivative nature of the higher-multipole vector-meson interactions. Following [40] we notice that the range of these delta contributions is shorter than that of proper vector meson exchange, from which we do not expect them to contribute significantly to the saturation of the C b , D b and D c couplings. This observation coincides with the findings of phenomenological studies of vector meson exchange in the charmed antimeson and charmed baryon systems -the pentaquarks -, which usually require to neglect the Dirac-delta contributions for the correct reproduction of the known pentaquark spectrum [48,49]. Ref. [48] points out how these contributions excessively distort the one pion exchange potential at relatively long distances, while Ref. [49] suggests that a possible rationale for neglecting the Dirac-deltas might be the exchange of heavier light mesons not explicitly included.
Going back to the relative importance of the components of the potential, it happens that W b , which will contribute to the D b coupling, is expected to be remarkably strong, where in Ref. [41] it was found to generate a large contribution to the binding energy of the Y (4230). Owing to the importance of the W b component to generate the Y (4230), at low momenta we will propose the hierarchy which might be representative of the expected hierarchy for the contact-range couplings D b and C a , C b , D c . It is worth noticing that the W c component has also been argued to receive a strong contribution from vector meson exchange, and that by calibrating its strength it is possible to reproduce the mass of the Y (4360) [41]. Thus, by repeating the previous arguments we arrive at which might indicate a strong D c coupling too. We warn that these hierarchies are merely qualitative in nature though, and that they are dependent on the assumption of which components of the potential generate binding for a given molecule. In particular, other orderings are possible. In a first approximation we will consider the first hierarchy we have proposed, i.e. the one given by Eq. (56) Finally it is important to stress that besides lightmeson exchanges, the EFT couplings will also receive contributions from the shorter-range cc components of the wave function of the Y states [47]. Owing to their shorter range, these contributions are naively expected to have a smaller impact on the value of the EFT couplings, which is why we do not discuss them directly. If we model the attraction coming from the cc component as charmonium exchange, as proposed in [27,[43][44][45], the size of their contributions to the contact-range couplings will scale as 1/m 2 J/ψ to be compared with 1/m 2 S and 1/m 2 V for scalar and vector mesons. Yet it is important to stress that within the EFT framework the specific origin of the contact-range couplings is inconsequential: they are determined from the same experimental information regardless of whether they come from light-meson exchanges or a compact charmonium core. From the molecular interpretation of the Y (4230) and the phenomenology we have previously discussed (derived from [41]), we propose the following power counting We will call it Y -counting, as it is constructed to reproduce the Y (4230) with a single parameter: D 0b . Indeed, in this counting the leading order (LO) potential reads We can determine the D 0b coupling from solving the bound state equation where φ(k) is the vertex function, and M mol the mass of the molecular state we are predicting. However, the contact-range potential is singular and has to be regularized (and renormalized, as explained in the next paragraph), which we do by including a separable regulator function where we will choose a Gaussian regulator is the unregularized LO potential, which is given by Eq. (59) within the Y -counting.
The LO potential reduces to a coupling constant that can be expressed as a linear combination of the subset of C Ia , C Ib , D Ib and D Ic couplings that are non-zero at LO. We denote this coupling as C LO (Λ) in the second line of Eq. (61). It happens that for the Gaussian regulator we have chosen, the bound state equation (Eq. (60)) can be solved analytically, yielding: with γ = 2µ(M th − M mol ) the wave number of the bound state, where the loop integral I(γ, Λ) reads with erfc (x) the complementary error function. In this analytic expression it can be appreciated that I(γ, Λ) depends on the sign of γ: for γ > 0 we will talk about bound states, which would correspond to the numerical solutions of Eq. (60) or to negative energy solutions of the Schrödinger equation with the asymptotic behavior Ψ( r) ∝ e −γr /r for r → ∞. In contrast, for γ < 0 we will have virtual states instead, i.e. negative energy solutions of the Schrödinger equation with the asymptotic behavior Ψ( r) ∝ e +γr /r for r → ∞. If γ → 0 both the bound and virtual state solutions lead to a two-body cross section of σ → 4π/γ 2 , which implies that virtual states (despite not being bound) might still be detectable provided they are close enough to threshold. For this reason we will also include virtual states in our predictions. For renormalizing the potential, the first step is to make the couplings in V LO C dependent on the cutoff Λ, i.e. we will have C LO = C LO (Λ) as already indicated in Eqs. (61) and (62). Alternatively, if we expand the potential in its components (Eq. (59)) what we will have is C Ia = C Ia (Λ), C Ib = C Ib (Λ), etc., that is, all couplings (or the particular subset that survives at LO in the power counting we are using, e.g. D 0b = D 0b (Λ) for the Y -counting) will depend on the cutoff. The second step is to determine the couplings from the condition of reproducing a series of observable quantities (for instance, the physical masses of a candidate molecular state) for a given cutoff Λ. From this procedure we expect the relative cutoff dependence of the observables predicted in the EFT to scale as O(Q/Λ) at LO. As a consequence by taking Λ ∼ M the cutoff dependence of the LO calculation will be of the order of the subleading order corrections. Alternatively, by varying the cutoff around the breakdown scale M it will be possible to estimate the uncertainty of the LO results, though in general this method is known to underestimate the errors and we will complement it with a second method for calculating the errors.
In our case, we determine (or renormalize) D 0b from by reproducing the Y (4230) mass in Eq. (60), for which we will use the central value in Eq. (1), i.e. M Y = 4218.6 MeV (the measurement of Ref. [19]). We will use the cutoff Λ = 0.75 GeV, which is around the rho meson mass. For estimating the uncertainties of our predictions, we will first allow the cutoff to move in the (0.5−1.0) GeV range, which are values very commonly used in the literature. This choice results in a coupling of for Λ = 0.75 (0.5 − 1.0) GeV. Yet, we will consider also the EFT uncertainties coming from the wave number of the Y (4230) (from our input), as well as the wave number of the bound states we are predicting (from our output). We will estimate the relative size of this error to approach γ/m ρ when γ is small, with γ = 2µ(M th − M mol ) the wave number of a given state and m ρ the rho meson mass. The final uncertainty in the two-body binding energy B of the states we predict will be given by where the first and second errors are the cutoff and in- and γ max = max(γ, γ Y ), with γ Y the wave number of the input state, i.e. the Y (4230). The specific form of the intrinsic error converges to the more common relative γ/m ρ error and its equivalent to it within the LO Yet, in situations where the convergence of the EFT is slow (as it is probably the case here) this choice is more suitable than the usual γ/m ρ ratio. The reason lies in the bound states for which the binding energy is close to the limit of convergence of the EFT: as the γ/m ρ ratio increases for a given state, the error will eventually be compatible with the state not binding, which curiously would have not happen if the state was less bound in the first place. Our choice avoids this problem, while acknowledging at the same time that there will be large uncertainties in the mass of such a state.
With D 0b we can now predict the full spectrum of charmed meson-antimeson molecules in the Y -counting. For this, we first determine the unregularized potential for each molecular configuration by adapting Table II to the Y -counting, where the resulting potentials are listed in Table III. We then regularize the potential with a Gaussian regulator and plug it into the bound state equation (Eq. (60)). For the masses of the S-and P-wave mesons, which are required for calculating M th and µ in Eq. (60), we have used the isospin average of the masses listed in the RPP [21]. The poles in the first (second) Riemann sheet correspond to the bound (virtual) states predictions, which we list in Table III. If we take into consideration the LO EFT potential in the Y -counting in more detail, we quickly realize that the 1 −− D sDs1 and 0 −+ , 1 −− and 2 −− D * sDs1 and 3 −− D * sD * s2 configurations are attractive. For the 1 −− D sDs1 system, we find a bound state at which is compatible with the mass and quantum numbers of the Y (4500) recently detected by BESIII [23]. For the Y (4230), we notice that the Y (4620) has been produced in e − e + [24,25] and thus we expect its quantum numbers to be J P C = 1 −− , where we find a pole with mass with the superscript V indicating a virtual state solution. This virtual state has the correct quantum numbers to be a plausible interpretation of the Y (4620). The problem is that virtual states are very difficult to observe unless they are close to threshold, and in this case the center-ofmass energy of this pole is 25 +43 −22 MeV below the D * sD s1 threshold. The uncertainties in the exact location of this virtual state are large though and it might happen that  it is closer to threshold than our first approximation in the Y -counting suggests. There are potentially important effects that we have not considered yet, including the role of the D 0c coupling or the coupled-channel dynamics between the D * sD s1 and D * sD * s2 channels, which are about 34 MeV away from each other. As we will see, these effects will change the pole into a virtual or bound state close to threshold.
Besides the 1 −− D * sDs1 system, it is also interesting to consider its 2 −− counterpart, as their potential is the same in the present counting, leading to the same mass prediction in the same Riemann sheet: However the most attractive configuration is the In Table III we predict the 3 −− D s2D * s system to be bound and have a mass of The prediction of three states (1 −− , 2 −− and 3 −− ) with a mass in the vicinity of (4620 − 4640) MeV will not change after the inclusion of the second non-diagonal contactrange interaction (D 0c ) and coupled-channel effects, as we will see in the following lines.
D. What about the Y (4360)?
Another charmonium-like state which has been suspected to be molecular is the Y (4360), with a mass of 4382.0 MeV according to its latest observation by BE-SIII [19]. However it is difficult to predict in the Ycounting, where the reason is that its LO potential reads which is too weak to form a deep enough bound state. There exists attraction though, with a virtual state being expected merely 2 MeV below threshold, check Table III. This might change if we promote a second coupling to LO. In view of the phenomenological model of Ref. [41] the D 0c coupling seems a good choice, with its promotion entailing a second power counting: which we will call the Y * -counting. In this counting the LO potential in the Y (4360) channel reads where the new coupling D 0c can be determined from the mass of the Y (4360). Repeating the same steps as with the Y (4230), we find for Λ = 0.75 GeV, where the values in parentheses correspond to the (0.5 − 1.0) GeV range.
The LO potential for the other meson-antimeson configurations and the predictions for this counting are shown in Table IV. It is interesting to notice the J P C = 1 −− D * s D * s1 channel, where the potential reads that is, half the strength of the potential in the Y (4360) channel. This time we predict a bound state at which is close to the experimental location of the Y (4620). Besides, there is a prediction of a J P C = 2 −− D * s D s1 virtual state close to threshold, with a mass and, as happened in the Y -counting, the most attractive hidden-strange molecule is, as usual, It is also worth noticing the existence of a 2 −+ D * s D * s2 virtual bound state with mass which one might be tempted to identify with the X(4630) state found by the LHCb [26], except that it is unlikely to be seen owing to its location in the complex plane. Being so deep as to be below the D * sD s1 threshold, this result only underscores that it might be advisable to include coupled channel effects, which we will do in the following lines.

E. Coupled channel effects
The mass difference of the ground and excited states of the P-wave charmed mesons is small in comparison with their S-wave counterparts. Indeed, we have that respectively, which happens to be comparable with the binding energies of the states we have been predicting in Tables III and IV. In terms of a momentum scale for the D * D * 2 and D * sD * s2 systems, the previous mass differences would correspond with a coupled-channel momentum of Λ CC = 288 and 275 MeV, which belongs to the light scales. This suggests that the inclusion of coupled channel effects would be a possible improvement to our EFT description.
The molecular configurations affected by coupled channel effects are the J = 1, 2 D * D 1 -D * D * 2 and D * sD s1 -D * sD * s2 systems. If we limit ourselves to the scope of the two power countings proposed in this work, the coupled channel potentials that we will obtain in these configurations for the non-strange sector are where the sign depends on the C-parity configuration, while for the hidden-strange sector the potentials will be half the previous ones The most important effect of these potentials is to favor the formation of bound states in the J = 1, 2 D * D 1 and D * sDs1 systems, while disfavoring their higher mass counterparts J = 1, 2 D * D * 2 and D * sD * s2 (which actually disappear). In Table V we summarize the new predictions for the D * D 1 and D * sD s1 systems when coupled channels are included in the Y -and Y * -countings.
The inclusion of these coupled channels does not modify the determination of the D 0b coupling in the Ycounting, only its predictions for the affected molecular configurations. The most significant changes happen in the 1 −− D * D 1 channel, where the new prediction is and the appearance of poles with positive C-parity, which might provide a plausible molecular interpretation for the X(4630) MeV .
The first configuration is a virtual state with respect to the D * s D s1 and D * s D * s2 thresholds, that is, the (II,II) Riemann sheet (denoted with the superscript V ′ ), from which we do not expect it to be observable in experiments. The second configuration generates a bound state basically at threshold, which might correspond with the X(4630), though depending on the cutoff it might become more or less bound, where in this later case it can move to the (II,I) Riemann sheet (i.e. becoming a resonance).
For the Y * -counting, the situation is slightly different. As a consequence of the new channel we have to recalibrate the D 0c coupling, resulting in   that is, a slightly weaker coupling (now the higher mass channel provides additional attraction to this system). In this counting the J P C = 1 −+ D * sD s1 configuration is more attractive than before, with the pole moving to the (I,II) Riemann sheet and located at where the V * superscript indicates the previous combination of sheets in the complex plane. Being in the (I,II) sheet, the 1 −+ state is still too deep to be observable and hence the Y * -counting with coupled channels still disfavors interpreting it as the X(4630). In contrast, the 2 −+ configuration still binds and happens to be located basically at threshold, again: where the errors are asymmetric as a consequence that for softer (harder) cutoffs it increases binding (becomes a virtual state below threshold). The bound state solution is not far away from the experimental mass of the X(4630) (and compatible once we take into account the experimental uncertainty, i.e. 4626 +18 −110 MeV), with its preferred spin-parity 1 − or 2 − . Thus the X(4630) can potentially correspond with our predicted 2 −+ D * sDs1 state.

F. Limitations of the present EFT approach
Both the Y -and Y * -counting are simplifications of the complex operator structure of the contact-range potential for this type of molecules. Despite having five candidates, only two of them have clear molecular interpretations, the Y (4230) and Y (4360); two of the hidden-strange candidates seemingly are the SU(3)-flavor partners of the previous two non-strange states, the Y (4500) and Y (4620) that might be partner states of the Y (4230) and Y (4360), respectively. Meanwhile, the hidden-strange X(4630) state is more difficult to interpret without coupled-channel dynamics.
The previous limitations, combined with the phenomenological observation that the E1 term (and to a lesser extent the M2 term) of the S-wave potential might be relatively strong, lead to the choices behind the Yand Y * -countings. These two countings are seemingly able to explain the five previous states, specially when coupled channel dynamics are included. Yet, they are incomplete. Luckily it is not difficult to analyze in which molecular configurations the previous two countings are more likely to fail.
The crucial assumption of the Y -and Y * -counting is the following ordering of the size of the couplings for which there is circumstantial evidence from saturating the value of these couplings with scalar and vector  meson exchange and then estimating the magnitude of the later from vector meson dominance or the condition of reproducing the Y (4230) and Y (4360) [41]. The obvious problem is that the previous ordering only works when the numerical factors in front of all the couplings is of O(1). Thus, the Y -and Y * -countings are expected to work well for the DD 1 , DD * 2 and J = 3 D * D * 2 systems, while their applicability in other configurations will be more limited.
This will generate systematic deviations from our predictions, which can be estimated up to a certain extent. For instance, the C Ia couplings are more likely to be attractive than repulsive and are also likely to be attractive in both I = 0, 1 configurations (with I = 0 being the most attractive one): The reasons behind this hypothesis are on the one hand the saturation of the couplings from scalar and vector meson exchange and on the other the fact that similar relations can be deduced from the molecular candidates in the D ( * )D( * ) system, e.g. the X(3872) and Z c (3900/4020), where caveats might apply as these are different systems and their couplings are not the same. If the previous relation happens to be the case, this indicates that we are underestimating the magnitude of the attraction in configurations where the coefficient in front of D 0b and D 0c is small and also in the hidden-strange configurations. Thus, this argument points towards more binding for the hidden-strange molecules and for most of the D * D 1 and D * D * 1 configurations. The next deviation regards the spin-spin coupling C 0b , for which the sign is not known for sure. There exist two hypotheses in the literature: which we will call scenarios A and B (in analogy with how these possibilities have been previously named in the pentaquark case [50]). Scenario A generates a spectrum in which the molecules with higher light-spin content are heavier, while scenario B predicts the opposite pattern. Scenario A is the standard one for compact hadrons, but it might also happen in isoscalar two-hadron states if the contact-range couplings are saturated by the direct interaction of the light-quarks within the hadrons, as proposed in Ref. [51]. Instead, if the contact-range couplings are saturated by light-meson exchange [40] (plus the assumption that the Dirac-delta terms appearing in the spin-spin vector-meson exchange terms can be neglected at the scales relevant for the formation of bound states, see e.g. the brief discussion after Eqs. (52)(53)(54)(55) or the ones in Refs. [48,49] for the pentaquark case), then the expected outcome will be scenario B. Which of these two scenarios is correct matters for the molecular interpretation of the X(4630), in particular whether it is a 1 −+ or 2 −+ state. The reason is that in these two channels the potential reads where for simplicity we have not specified the isospin or flavor representation. If scenario A is correct, the 1 −+ configuration (which is already attractive, only not enough as to generate a pole near threshold) will be more attractive than expected, while the contrary will happen for the 2 −+ configuration. This might potentially change our conclusions and indicate that the most promising molecular interpretation of the X(4630) is a 1 −+ D * sDs1 bound state, which will be in agreement with Ref. [28]. Conversely, if scenario B is correct, this will further cement the conclusion that the X(4630) is better reproduced as a 2 −+ D * s D s1 bound state. In the next lines we will explicitly consider two phenomenological models as a way to further understand the molecular spectrum and its relation with the factors we have discussed here.

G. Comparison with phenomenology
Actually, there is a possible way to address the unknown factors that we have discussed in the previous subsection: to calculate the molecular spectrum in a phenomenological model and compare these predictions with the EFT ones.
Even though phenomenology is intrinsically modeldependent, it nonetheless contains all the interaction structures we have discussed (V a , V c , W b , W c or their contact-range EFT counterparts C a , C b , D b , D c ), including those that we are currently unable to determine owing to the absence of suitable molecular candidates (i.e. C a and C b ). In this regard, if EFT and phenomenological predictions happen to converge, that might indicate that we are on the right track for determining the spectrum of S-and P-wave charmed meson molecules.
For this comparison we have chosen the following two phenomenological models: (i) The one boson exchange (OBE) model [52,53], in which the potential among hadrons is the sum of the potential generated by the exchange of a few light-mesons, in our case the scalar and vector meson exchange potential of Eqs. (52)(53)(54)(55), but modified by the inclusion of finite-size effects (i.e. form factors).
(ii) A saturation model in which the finite-range potential generated by the exchange of a light meson is simplified to a contact-range potential with a coupling strength proportional to that of the original finite-range potential. By means of a renormalization group (RG) equation, the previous coupling strength is in turn modified as to take into account the different ranges of the scalar and vector mesons [41], which is why we will refer to this model as the RG-saturation model.
We describe the technical details of these models in Appendix C, while here we concentrate on comparing their results to those of EFT. Yet, we will briefly comment on a few commonalities of the two models. First, in the OBE model there is an unknown form factor and cutoff, while in the RG-saturation model there is an unknown proportionality constant for the saturated coupling. In both cases we determine these unknown factors by reproducing the X(3872) as a J P C = 1 ++ , I = 0 D * D molecule (that is, we are assuming a molecular X(3872)). Second, while most of the couplings of the scalar and vector mesons are relatively well known, this is not the case for the E1 and M 2 components of the vector meson potential, i.e. the f ′ V and h ′ V couplings of Eqs. (54) and (55). Here we fix f ′ V and h ′ V from the condition of reproducing the masses of the Y (4230) and Y (4360).
We also provide error estimations for the OBE and RG-saturation models. The status and parameters of the scalar meson are much more uncertain than those of the vector mesons. In particular, the scalar meson mass is not well known, where the RPP cites the m S = (400 − 550) MeV range [21]. This uncertainty, which by itself implies a variability of the overall strength of the potential of the order of (30 − 40)%, can be propagated to the predictions of both models. For the RGsaturation model a second source of uncertainty is that it uses the contact-range approximation, which is a good approximation only if the predicted bound state is sufficiently large in size. Otherwise the wave function of the bound state will be able to probe the short-range details of the interaction binding them. For this we will consider the same type of γ/m V relative error that we already included in the EFT calculations, i.e. [(1 + γmax mV ) ±1 − 1] with γ max the largest of the relevant binding momenta (either of the input states, i.e. the Y (4230), or the predicted states), check Eqs. (65) and (66) and the related discussion for details.
The comparison among these two models and the two power countings we have discussed is shown in Tables VI and VII for the non-strange and hidden strange sectors. It is important to notice that all calculations include the D 1D * -D * 2D * and D s1D * s -D * s2D * s coupled channel effects when relevant. From inspecting Tables VI and VII, it can be appreciated that the two phenomenological models tend to predict a very similar set of molecular states. There are exceptions though, the most notable one being the J P C = 1 −+ D 1D * and D s1D * s configurations, which are predicted as unbound in the OBE model but as bound in the RG-saturation model. This configuration is interesting as it might correspond to the X(4630), for which the spin is either J = 1 or 2, and in this regard the RG-saturation model predicts a candidate state for both spins.
For the rest of the states, the predictions of both phenomenological models basically coincide with the ones of the Y * -counting, while there are disagreements with the Y -counting for the D * 2D and D * s2Ds configurations, for the simple reason that the coupling providing binding for these molecules is missing, and also the J = 0 D 1D * and D s1D * s molecules. In this last case, the Yand Y * -countings predict binding for the J P C = 0 −+ and 0 −− , respectively. However, as we explained in the previous subsection, the predictions for these configurations are uncertain, particularly in the Y -counting. The comparison with the two phenomenological models presented here supports the previous conclusion, except that for the Y * -counting the predictions basically align with phenomenology, and also with the recent prediction of a J P C = 0 −− D 1D * state in Ref. [9]. Yet, caution is advised in view of the discrepancies in the prediction of its mass in the OBE and RG-saturation models and EFT, suggesting large uncertainties.

V. CONCLUSIONS
In this manuscript we have considered the Y (4230), Y (4360), Y (4500) and Y (4620) resonances from a molecular perspective, where the interaction between the ? ? charmed meson and antimeson is constrained by their heavy-and light-quark symmetries, i.e. HQSS and SU(3)-flavor. Though it is unlikely that the Y (4230) or any other of the aforementioned resonances is a purely molecular state, it is also not particularly probable that they are completely non-molecular. For simplicity we have treated these states as pure molecules, without admixture from shorter-range charmonium components, and then built an EFT description of the potential between their hadronic components, which at leading order only involves contact-range interactions.
The resulting EFT potential contains four unknown couplings per isospin or SU(3)-flavor representation, i.e. eight couplings in total. This implies that there are not enough molecular candidates available to determine these couplings, a limitation that we propose to overcome by invoking a series of simplifications. On the basis of a preexisting light-meson exchange phenomenological model [41], we are motivated to make the assumption that the interaction between the S-and P-wave charmed mesons conforming the Y (4230) state is dominated by a non-diagonal DD 1 → D 1D term involving an electric-like dipolar E1 D → D 1 transition mediated by vector-meson exchange. In this case, which we call the Y -counting, a molecular Y (4230) will imply the existence of a few partners, among which there are states with mass and quantum numbers potentially compatible with those of the Y (4360), Y (4500) and Y (4620). The Y (4360) can be interpreted as a J P C = 1 −− D * D 1 bound state, the Y (4500) as a J P C = 1 −− D sDs1 one and the Y (4620) as a J P C = 1 −− D * sDs1 virtual state, with this later state probably having 2 −− D * sD s1 and 3 −− D * sD * s2 partners of similar mass.
The previous predictions can be further refined by considering a second non-diagonal D * D 1 → D 1D * term involving a magnetic-like quadrupolar M2 contribution to the D * → D 1 transition. We call this possibility the Y *counting, as we determine the new coupling associated with the M2 interaction from the mass of the Y (4360). The inclusion of this interaction improves the predictions of the mass of the candidate molecular configurations for the Y (4620) -the 1 −− D sDs1 configuration -which now becomes a bound state.
Yet, a further improvement is possible: the mass difference of the ground and excited states of the P-wave charmed mesons is merely 39 MeV (34 MeV) for the nonstrange (strange) case, suggesting the inclusion of the coupled channel effects related to the D 1 → D * 2 transitions. This choice results for instance in a moderate improvement of the attraction and binding of the candidate configuration for the Y (4620). The molecules where coupled channel effects are more conspicuous are the J P C = 1 −+ and 2 −+ D * sD s1 systems, which are both attractive now. In particular, there is a near-threshold 2 −+ D * sD s1 bound state that might very well correspond to the X(4630) discovered by the LHCb [26], the spin and parity of which are expected to be J P = 1 − or 2 − . However there are limitations to the conclusions we draw in this manuscript. For starters neither the Y (4230) or Y (4360) are expected to be fully molecular, particularly when one considers their relatively large binding energies, which suggests that the cc components of their wave functions [11,47] are probably important. Yet, we believe that the symmetry considerations on which our analysis rests will more likely than not supersede the specific details about the structure of these two resonances: the Y (4500) and Y (4620) are actually predicted as hidden-strange partners of the Y (4230) or Y (4360). This point could be elucidated by future theoretical analyses.
There is a second caveat regarding the two power countings built for deriving our conclusions, which are based on the assumption that the non-diagonal E1 and M2 interactions are considerably stronger than the central and spin-spin ones. Even though there might be reasons to think this is the case, the numerical factors in front of these E1 and M2 terms change considerably for the different charmed meson-antimeson configurations, particularly when the total spin is low. As a consequence these power countings are expected to work better in molecules such as DD 1 , DD * 2 and J = 3 D * D * 2 (for which the numerical factor in front of the E1 and M2 couplings -D 0b and D 0c in this work -is close to one), but not so well for the other configurations of these systems. This means that it is more likely than not that we are underpredicting the amount of attraction of a few configurations of the D * D 1 and D * sD s1 systems. Of particular importance here is the sign of the spin-spin coupling: an attractive spin-spin term will reinforce the conclusion that the X(4630) is a 2 −+ D * sDs1 molecule, while a repulsive one will favor binding in the 1 −+ D * sD s1 configuration, which will thus become the preferred interpretation of the X(4630) (as happens in [28,30]). Here we have calculated the molecular spectrum in two phenomenological models, both of which support the idea of binding in the 2 −+ D * sDs1 system, but also showing that the 1 −+ D * sD s1 molecule could bind as well. The eventual experimental observation of molecular candidates in these systems together with improved lattice calculations (in the line of Refs. [54,55]) or more refined phenomenological models might be used to overcome the ambiguity of the previous results in the future.

ACKNOWLEDGMENTS
We would like to thank Xiang-Kun Dong  ( * ) 1(2)D ( * ) molecules. In particular we are interested in whether OPE is perturbative or not: previous experience in the D ( * )D( * ) two-body system, where OPE is indeed perturbative [32,34], suggests that this is likely to be the case too when one of the S-wave charmed mesons is substituted by a P-wave charmed meson. Yet, owing to the differences between these two systems, this hypothesis should be explicitly checked. As we will see, OPE is most probably still perturbative when P-wave charmed mesons are involved.
First, we begin by noticing that there are two types of OPE for the D ( * ) 1(2)D ( * ) molecules: diagonal (e.g D * D * 2 → D * D * 2 ) and non-diagonal (e.g D * D * 2 → D * 2D * ). For the diagonal one, the relevant Lagrangians in the superfield notation are which are basically the non-relativistic limit of the Lagrangians in [42], and where g 1 and g ′ 1 are the axial coupling constants, H and T k the S-and P-wave charmed meson superfields and A represents the non-relativistic axial pion current, which we define as with f π ≃ 132 MeV the weak pion decay constant, τ a the Pauli matrices in isospin space, π a the pion field and a = 1, 2, 3 an isospin index. If we choose the subfield notation instead, then the Lagrangians will read with q s , q p , σ L and S L having the same meaning as in Eq. (25). There is also the non-diagonal pion interaction, which reads in the standard superfield notation [42], with Λ χ a typical chiral energy scale (usually chosen to be around 1 GeV) and where the difference between the terms proportional to the h 1 and h 2 couplings lies in how the subindices of the T superfield are contracted, i.e. T · ∂ and T · A for h 1 and h 2 , respectively. If we write this Lagrangian in the subfield notation, we will obtain instead where the second term (proportional to the coupling f ) actually cancels out (because ∂ i A j is proportional to ∂ i ∂ j π a and thus symmetric). The relation of the couplings h and f with those of Eq. (A6) is With the previous Lagrangians we can derive the OPE potential in momentum space, leading to where q refers to the exchanged momentum, m π ≃ 138 MeV the pion mass and µ π the effective pion mass for the non-diagonal transitions (µ 2 π = m 2 π − ω 2 π , with ω π = m(D * 1(2) ) − m(D ( * ) )). This potential is not purely S-wave, but contains tensor forces with angular momenta of L = 2 for the diagonal piece and L = 2, 4 for the nondiagonal one. These tensor forces require S-to-D-and Sto-G-wave mixing to operate. Here we will assume that the higher partial wave components of the molecular wave functions are kinematically suppressed and consequently ignore these tensor forces. The S-wave, finite-range piece of the OPE potential can be easily isolated and reads where the dots indicate contact-range contributions or tensor forces. For determining the momentum scale above which OPE becomes non-perturbative for S-wave scattering, which we will call Λ S , we will rewrite the previous contributions as where µ is the reduced mass of the system,Ô L12 refers to the specific spin-spin operator that applies and Λ OPE is the characteristic OPE momentum scale (which is proportional to Λ S , the scale we are looking for). In this form, the ratio between iterated and tree-level OPE is given by where the operators have been evaluated for an S-wave scattering state, V S is the S-wave projection of the OPE potential, G 0 (E) = 1/(E − H 0 ) is the resolvent operator (with E is the center-of-mass energy and H 0 the free the particular two charmed meson state under consideration. In this second case we obtain for J = 1 , Finally, we consider the quadrupolar operators and their scalar products, the evaluation of which yields (in terms of light spin) from which, as previously mentioned for Σ † L · Σ L , one can derive the matrix elements in the charmed meson basis provided we have the heavy-and light-spin decomposition of the two-body states. If we use instead Eqs. (36)(37)(38)(39), we can directly write the matrix elements of the projection of the quadrupolar operators in the charmed meson basis as follows from which we will be able to derive the numerical factors found in Table II for the D c coupling.
Finally, we consider coupled-channel effects for the J = 1, 2 molecules in the bases where the correct C-parity combinations are implicitly assumed (though not written down). In the previous bases, the matrix elements of the non-diagonal operators all of which change sign if the C-parity is changed. For completeness, though it is not strictly a transition operator, it is useful to provide too the matrix elements of the spin-spin operator in the previous two bases: and vector meson exchanges in Eqs. (52)(53)(54)(55), yet within the OBE model these potentials are modified as a consequence of the finite size of the hadrons involved. For this, a form factor and a cutoff are included as follows: V M(finite) ( q , Λ OBE ) = V M(point−like) ( q ) f 2 M ( q) , (C1) with V M(point−like) the potential for point-like hadrons, i.e. Eqs. (52)(53)(54)(55), V M(finite) the potential for finite-size hadrons and where f M (x) is the form factor representing the finite-size effects for the exchange of meson M . In particular we will consider multipolar form factors of the type where m is the mass of the exchanged meson, q 2 = −q 2 0 + q 2 the square of the 4-momentum, Λ M a cutoff and n F and exponent. In principle it is possible to use a different form factor and cutoff for every of the exchanged mesons, which is what is done when the OBE is applied in the twonucleon sector [52,53]. But here we will limit ourselves to a unique form factor and cutoff (otherwise there will be not enough input data to determine the cutoffs).
Once the cutoff and the exponent have been chosen predictions are easy to make by solving the resulting potential in the Schrödinger equation. For the exponent we will choose n F = 3, while the cutoff is chosen by applying the OBE model to the X(3872) as a 1 ++ D * D molecule and then determining the cutoff from the condition of reproducing the mass of the X(3872). This will lead to Λ OBE = 1.56 GeV for the parameters we specify below.
Besides the form-factor and the cutoff, the parameters of the OBE model are the masses and couplings of the light mesons that are exchanged. For the masses we will use m σ = 475 MeV, m ρ = 770MeV and m ω = 780 MeV. For the couplings we begin with the ones appearing in the diagonal components of the potential, Eqs. (52) and (53), for which we take g σ = 3.4 (derived from the linear sigma model [64] and the quark model [65]), g V = 2.9 (from the universality of the vector meson coupling [66][67][68], which requires g V = m V / √ 2f π with f π = 132 MeV), f V = g V κ V i with κ(D ( * ) ) = 2.9 and κ(D ( * ) 1(2) ) = 2κ(D ( * ) ) [41]. For the non-diagonal components of the potential, Eqs. (54) and (55), we define f ′ V = g V κ E1 and h ′ V = g V κ M2 and determine κ E1 and κ M2 from the condition of reproducing the masses of the Y (4230) and Y (4360) as 1 −− DD 1 and D * D 1 molecules. This results in κ E1 ≃ 4.9 and κ M2 ≃ 22.8.

The RG-improved saturation model
The second phenomenological model we will consider is that of a contact-range theory in which the strength of the couplings constants is derived from the condition that they are saturated by light-meson exchanges [69,70]. The idea is that at the low momenta characteristic of a hadronic bound state the details of the light-meson exchange potential are not resolved and thus its effects can be reduced to a contact-range interaction of suitable strength. For a generic light-meson exchange potential of the type where the subindex M indicates the particular meson considered, v M is a prefactor containing coupling constants and operators and the dots denote either Diracdeltas or higher angular momentum components, the equivalent saturated coupling will be proportional to where Λ sat refers to the regularization scale used to regularize the contact-range potential, which is expected to be of the order of the mass of the exchanged light-meson. The proportionality constant is in principle unknown. It is also important to stress that the previous refers to the specific saturation method proposed in [40], which explicitly ignored the contact-range pieces of the light-meson exchange potentials -the dots in Eq. (C3) -as they represent physics at scales corresponding to the hadron internal structure (instead of the light-meson exchange range, which is what interest us). When considering the saturation from the scalar and vector mesons, there is the problem of how to combine them properly: the mass of the scalar and vector mesons, though of the same order of magnitude, are different enough as to justify different values of Λ sat . What we will do instead is to combine their contributions by means of the RG equations proposed in [71], which in practical terms implies the following rule for combining the saturation from the scalar and vector mesons with m S and m V the scalar and vector meson masses, respectively, and α an exponent related with the behavior of the wave function at the distances m V r = O(1) and m S r = O(1). For a detailed explanation we refer to [41], where it was also proposed the value α = 1 for the exponent. With all this, the saturated contact-range coupling for the D ( * )D ( * ) 1(2) molecules is given (modulo a proportionality constant) by where κ V = f V /g V , κ E1 = f ′ V /g V , κ M2 = h ′ V /g V and with m S , m V , µ V , ω V and M as defined in Eqs. (52)(53)(54)(55). The proportionality constant can be determined from the condition of reproducing the mass of a given molecular state. For this, the saturated contact-range coupling is to be regularized as usual where Λ sat ∼ m V is the saturation cutoff and f (x) a regulator function. Then we plug this potential into the bound state equation, Eq. (60), in order to determine the masses of the molecular states. For concreteness we take Λ sat = 1.0 GeV, f (x) = e −x 2 and use the X(3872) as a 1 ++ D * D bound state as the reference state from which to fix the proportionality constant. For the masses and couplings we use the same values as for the OBE model except for κ E1 and κ M2 , which we determine again from reproducing Y (4230) and Y (4360) as 1 −− DD 1 and D * D 1 bound states, resulting in κ E1 ≃ 3.9 and κ M2 ≃ 24.5.