X(3872) is not a true molecule

A solvable coordinate-space model is employed to study the $c\bar{c}$ component of the X(3872) wave function, by coupling a confined $^3P_1$ $c\bar{c}$ state to the almost unbound $S$-wave $D^0\bar{D}^{*0}$ channel via the $^3P_0$ mechanism. The two-component wave function is calculated for different values of the binding energy and the transition radius $a$, always resulting in a significant $c\bar{c}$ component. However, the long tail of the $D^0\bar{D}^{*0}$ wave function, in the case of small binding, strongly limits the $c\bar{c}$ probability, which roughly lies in the range 7-11%, for the average experimental binding energy of 0.16 MeV and $a$ between 2 and 3 GeV$^{-1}$. Furthermore, a reasonable value of 7.8 fm is obtained for the X(3872) r.m.s. radius at the latter binding energy, as well as an $S$-wave $D^0\bar{D}^{*0}$ scattering length of 11.6 fm. Finally, the $\mathcal{S}$-matrix pole trajectories as a function of coupling constant show that X(3872) can be generated either as a dynamical pole or as one connected to the bare $c\bar{c}$ confinement spectrum, depending on details of the model. From these results we conclude that X(3872) is not a genuine meson-meson molecule, nor actually any other mesonic system with non-exotic quantum numbers, due to inevitable mixing with the corresponding quark-antiquark states.


Introduction
The Belle Collaboration discovered [1] the charmoniumlike state X(3872) almost a decade ago, observing it in the decay B ± → K ± π + π − J/ψ, with a significance in excess of 10σ. The new meson was then confirmed by CDF [2], D0 [3], BABAR [4], and recently also by LHCb [5]. In the 2012 PDG tables [6], X(3872) is listed with a mass of 3871.68 ± 0.17 MeV and a width Γ < 1.2 MeV. As for the quantum numbers, the possibilities are J P C = 1 ++ and J P C = 2 −+ according to the PDG [6]. The positive C-parity assignment resulted from an analysis of angular distributions by CDF [7]. Recently, observation of the decay to γJ/ψ by Belle [8] unambiguously confirmed C = +. For further experiments on X(3872) production, see Ref. [6].
On the theory side, the discussion about the nature of X(3872) continues most vivid. Although the PDG [6] lists 1 ++ and 2 −+ as the meson's possible quantum numbers, and a recent BABAR analysis [9] of the 3π invariant mass distribution in the decay X(3872) → ωJ/ψ even seems to slightly favor the 2 −+ assignment, most model builders describe the state as an axial vector. For instance, model calculations of semi-inclusive B → η c2 + X processes [10] as well as electromagnetic η c2 decays [11] have been shown to disfavor the 2 −+ scenario. The same conclusion was reached in a tetraquark description of X(3872) [12], while pion exchange in a molecular picture would be repulsive in this case [13] and so inhibitive of a bound state. Finally, unquenching a 1 1 D 2 cc state by including meson-meson loops can only further lower the bare mass, which lies in the range 3.79-3.84 GeV for all quenched quark models we know of, thus making a 2 −+ charmonium resonance at 3.872 GeV very unlikely [14]. For further information and more references concerning X(3872), see e.g. a recent review [15], as well as our prior coupled-channel analysis [14].
The first suggestion of possible meson-meson molecules bound by pion exchange, in particular a DD * 1 state with quantum numbers 1 ++ or 0 −+ , was due to Törnqvist [13]. With the discovery of X(3872) just below the D 0 D * 0 threshold, this idea was revived, of course. In the present paper, we intend to study the issue, not from Törnqvist's pion-exchange point of view, but rather as regards its possible implications for models based on quark degrees of freedom. In this context, it is worthwhile to quote from Ref. [16], in which a molecular interpretation is advocated (also see Ref. [17]): "Independent of the original mechanism for the resonance, the strong coupling transforms the resonance into a bound state just below the two-particle threshold if a > 0 or into a virtual state just above the two-particle threshold if a < 0. If a > 0, the bound state has a molecular structure, with the particles having a large mean separation of order a." (Note that, here, a is the S-wave scattering length.) Also: "In this case [1 ++ ], the measured mass M X implies unambiguously that X must be either a charm meson molecule or a virtual state of charm mesons." In face of these peremptory claims about the molecular picture, we believe it is of utmost importance to study in detail the X(3872) wave function for a model in which the mechanism generating the meson is quark confinement combined with strong decay. Thus, we employ a simplified, coordinate-space version of the coupled-channel model used recently [14] to describe X(3872) as a unitarized and mass-shifted 2 3 P 1 charmonium state. The model's exact solvability then allows to obtain analytic expressions for the wave-function components, and follow bound-state as well as resonance poles on different Riemann sheets.
The model is outlined in Sec. 2, with details moved to Appendices A and B. Section 3 is devoted to S-matrix pole trajectories as a function of the two free parameters. In Sec. 4 the two-component wave function is analyzed for several parameter sets, and is then used in Sec. 5 to compute cc probabilities and root-mean-square (r.m.s.) radii. The dynamical vs. confinement nature of the poles is discussed in Sec. 6. Conclusions are drawn in Sec. 7.
2 The coupled cc -D  D *  system In our previous work on X(3872) [14], we described the state as a unitarized radially excited 1 ++ charmonium resonance, by employing the resonance-spectrum-expansion (RSE) [18] formalism. The RSE description of X(3872) amounted to coupling the most relevant meson-meson channels allowed by the Okubo-Zweig-Izuka (OZI) rule to a spectrum of bare 1 ++ cc states, and also the OZI-suppressed but experimentally observed [6] ρ 0 J/ψ and ωJ/ψ channels. Thus, we found that unquenching shifts the bare 2 3 P 1 state more than 100 MeV downwards, which settles as a very narrow resonance slightly below or on top of the D 0 D * 0 threshold. However, the RSE approach does not allow to obtain wave functions in a straightforward fashion. So for the purpose of the present paper, we resort to the equivalent [19] coordinate-space coupled-channel formalism of Ref. [20], which was used to study the influence of strong decay channels on hadronic spectra and wave functions, besides several more specific phenomenological applications. Furthermore, since here we do not aim to describe all aspects of the X(3872) resonance in a detailed way, but rather want to focus on the importance of the charmonium component in the wave function, we restrict ourselves henceforth to a simple two-channel system, viz. a 3 P 1 cc state coupled to an S-wave D 0 D * 0 channel, the latter being the dominant one, observed in the D 0 D 0 π 0 mode [6]. A partial study of this problem was already carried out in Ref. [21]. Consider now a system composed of a confined qq channel coupled to a meson-meson channel M 1 M 2 . Confinement is described through a harmonic-oscillator (HO) potential with constant frequency ω, having spectrum where ν is the radial quantum number, l c the qq orbital angular momentum, and m q = mq = 2µ q the constituent quark mass. In the scattering channel, no direct interactions between the two mesons are considered, with µ f and l f being the reduced two-meson mass and orbital angular momenta in the free channel, respectively. Transitions between the two channels are modeled via an off-diagonal delta-shell potential with strength g, which mimics string breaking at a well-defined distance a. The corresponding Hamiltonian, transition potential, and 2 × 2 matrix Schrödinger equation are given in Eqs.
(2)-(5), with the usual definition u(r) = rR(r), where R(r) is the radial wave function. The exact solution to these equations is derived in Appendix A: Once the 1 × 1 S matrix (cf. Eq. (19)) has been constructed from the wave function, possible bound or virtual states as well as resonances can be searched for. These occur at real or complex energies for which S blows up, or equivalently, when cot δ l f (E) = i (cf. Eq. (18)). Thus, unlike the pure HO spectrum, which is a real and discrete set of energies, the "unquenched" spectrum given by Eq. (5) includes complex energies, too, some of which even have no obvious connection to the original, "bare" levels. The corresponding poles, called dynamical, are the result of attraction in the meson-meson channel generated by the interaction with intermediate quark-antiquark states. The light scalar mesons f 0 (500), K * 0 (800), f 0 (980), and a 0 (980) [6] are archetypes of such resonances [22]. On the other hand, we designate by confinement poles the ones that can be linked straightforwardly to the bare states. Nevertheless, we shall show below -as we already demonstrated in previous work -that these two types of resonances are not so distinct after all, since minor parameter changes may transform one kind into the other. But regardless of the type of pole, the corresponding radial quantum number in Eq. (6) is related to the energy of the coupled-channel system by Only in the uncoupled case, that is, for g = 0, one recovers the original HO spectrum, with ν = 0, 1, 2, .... Now we apply the formalism to the coupled cc-D 0 D * 0 system. The cc channel is assumed to be in a 3 P 1 state, i.e., with l c = 1, implying the D 0 D * 0 channel to have l f = 0 or 2. Nevertheless, we shall restrict ourselves here to the S-wave channel only, which will be strongly dominant, especially near threshold. The fixed parameters are given in Table 1, where the meson masses are from the PDG [6], while the HO frequency ω and the constituent charm quark mass m c are as originally determined in Ref. [23] and left unaltered ever since. Thus, from Eq. (1) we get the lowest two HO states at E 0 = 3599 MeV and E 1 = 3979 MeV, respectively. The former should give rise -after unquenching -to the 1 3 P 1 charmonium state χ c1 (1P ) [6], with mass 3511 MeV, while the latter is the bare 2 3 P 1 state, which cannot so easily be linked to resonances in the PDG tables, though both X(3940) and X(3872) are possible candidates, in view of their mass and dominant DD * decay mode [6]. However, X(3940) may just as well be the, so far unconfirmed, 2 1 P 1 (1 +− ) state h c (2P ) [14].
The two remaining parameters, viz. the string-breaking distance a and the global coupling g, have to be adjusted to the experimental data. Nevertheless, these parameters are not completely free, as they both have a clear physical interpretation, albeit of an empirical nature. Thus, a is the average interquark separation at which 3 P 1 quark-pair creation/annihilation is supposed to take place, while g is the overall coupling strength for such processes. Note that we do not assume a particular microscopic model for string breaking inspired by QCD, like e.g. in a very recent paper [24]. Still, the values of a found in the present work are in rough agreement with our prior model findings, and even compatible [25] with a lattice study of string breaking in QCD [26]. Concretely, we have been obtaining values of a in the range 1-4 GeV −1 (0.2-0.8 fm), logically dependent on quark flavor, since the string-breaking distance will scale with the meson's size, being smallest for bottomonium. As for the coupling parameter g, its empirical value will depend on a, but also on the set of included decay channels. In realistic calculations, values of the order of 3 have been obtained (see e.g. our previous paper [14] on X(3872)).

Poles
The crucial test the present model must pass is its capability of generating a pole near the D 0 D * 0 threshold. Indeed, a dynamical pole is found slightly below threshold for different combinations of the free parameters a and g, several of which are listed in Table 2. Examples are here given of bound states, virtual bound states, and below-threshold resonances, the latter ones only occurring for S-wave thresholds as in our case. Note that poles of both virtual bound states and resonances lie on the second Table 1. Fixed parameters. Riemann sheet, i.e., the relative momentum has a negative imaginary part. From this table we also observe that larger and larger couplings are needed to generate a pole close to threshold when a approaches the value 3.5 GeV −1 .
We shall see below that this is due to the nodal structure of the bound-state wave function. Although a dynamical pole shows up near the D 0 D * 0 threshold, there still should be a confinement pole connected to the first radial 3 P 1 excitation at 3979 MeV. Well, we do find such a pole, for each entry in Table 2. In Table 3 a few cases are collected, with the parameters tuned to generate a dynamical pole at precisely the X(3872) PDG [6] mass of 3871.68 MeV. Note, however, that the associated confinement pole is not necessarily of physical relevance, since at the corresponding energy several other strong decay channels are open, which no doubt will have a very considerable influence and possibly even change the nature of both poles. As a matter of fact, in our prior paper [14], with all relevant two-meson channels included, the X(3872) resonance was found as a confinement pole, whereas dynamical poles were only encountered very deep in the complex energy plane, without any observable effect at real energies. So here we show these results only to illustrate that pole doubling may occur when strongly coupling S-wave thresholds are involved, as we have observed in the past in the case of e.g. the light scalar mesons [22] and D * s0 (2317) [27]. The issue of confinement vs. dynamical poles will be further studied in Sec. 6.
In order to better understand the dynamics of the different poles, we plot in Fig. 1 pole trajectories in the complex energy plane as a function of the coupling constant g, and for three different values of a. For vanishing g, the dynamical pole acquires a negative infinite imaginary part and so disappears in the continuum, whereas the confinement pole moves to the real energy level of the bare 2 3 P 1 state, i.e., 3979 MeV. As g increases, and for both a = 2.0 GeV −1 and a = 3.0 GeV −1 , the dynamical pole moves    Table. 3 are here marked by * ; (ii) arrows along curves indicate increasing g.  to the real axis below threshold, becoming first a virtual bound state and then a genuine bound state. Note that, in the latter case, the real part twice attains the X(3872) mass even before the pole reaches the real axis, but the corresponding imaginary parts are much too large as compared with experiment [6], so only the bound state can be considered physical. Finally, for a = 3.5 GeV −1 the pole does never reach the real axis, which would require an infinite coupling. For the other parameter sets listed in Table 2, we find intermediate situations. Another feature we can observe for all trajectories is an initial attraction and subsequent repulsion between the dynamical and the confinement poles.

Wave function
Now we are in a position to study the X(3872) boundstate wave function in several situations. We choose two values for the string-breaking parameter, viz. a = 2.0 GeV −1 and a = 3.0 GeV −1 . In Table 4 five different binding energies (BEs) are chosen with respect to the D 0 D * 0 channel, including the PDG [6] value labeled by X. We have computed and normalized (see Appendix A) the two-component radial wave function R(r) for each of the five cases. In Fig. 2 we depict the cases labeled by A, X and D, the other two representing intermediate situations. General features we immediately observe are the typical S-wave behavior of the D 0 D * 0 wave-function component R f , while the cc wave function R c is in a P state, the latter also having a node, as it is dominantly a first radial excitation. Furthermore, |R f | is larger than |R c | in most situations, for all r, except for unphysically large BEs (cf. plot D). Nevertheless, the two components are of comparable size for intermediate r values. Then, as the BE becomes smaller, the tail of R f grows longer, as expected, whereas R c always becomes negligible for distances larger than roughly 11-12 GeV −1 . Now, the increased R f tail affects the normalization of both R c and R f . Thus, the ratio |R f (r)|/|R c (r)| is quite robust for most r values, as it does not significantly change with the BE.
5 Probabilities and r.m.s. radii Having derived the X(3872) wave function for several scenarios, we can now straightforwardly compute the relative probabilities of the cc and D 0 D * 0 components (see Appendix A), with the results given in Table 5, for the five BEs and two a values from Table 4. Note that the probability in the D 0 D * 0 channel is only computed for normalization purposes, since in a more realistic calculation at least the D ± D * ∓ component would acquire a nonnegligible probability as well, as the corresponding threshold lies only 8 MeV higher. Nevertheless, our simplification is unlikely to have an appreciable effect on the cc probability and will only increase that of the D 0 D * 0 component accordingly. Also note that the cc probability includes all 3 P 1 states, with the 2 3 P 1 being dominant, because the corresponding bare eigenstate lies only 100 MeV higher. However, also the 1 3 P 1 state is nonnegligible in the physical X(3872) wave function. In the coupled-channel approach of Ref. [28], a 1 3 P 1 admixture of about 15% was found. Notice that -inevitably -unquenching not only mixes meson-meson components into the total bound-state wave function, but also quark-antiquark components of confinement states other than the one under consideration (also see Ref. [20]). Here, for a BE of 0.16 MeV, corresponding to the physical [6] X(3872), case X in Table 5, has a 7.48% cc probability for a = 2.0 GeV −1 and 11.18% for Fig. 2. Normalized two-component radial wave function R(r) for three BEs, corresponding to labels A, X, D in Table 4, and two a values. Upper curves: R f (r); lower curves: Rc(r). Left: a = 2 GeV −1 ; right: a = 3 GeV −1 .  Table 6. R.m.s. radii of the wave function, expressed in fm, for the cases specified in Table 4. This would then correspond to a cc probability roughly midway in the range 7.48%-16.98% (a = 2.0 GeV −1 ) or 11.18%-24.65% (a = 3.0 GeV −1 ). In the limiting case of zero binding, the cc probability would eventually vanish. Also notice that, in all five cases of Table 5, the cc probability rises by about 50% when a is increased from 2.0 to 3.0 GeV −1 . Nevertheless, if we take a = 2.0 GeV −1 as in our Ref. [14], we get a cc probability of 7.48%, very close the 7% found in Refs. [29,30]. Next we use the normalized wave functions and Eq. (17) to compute the X(3872) r.m.s. radius for the five cases discussed before (see Table 4), with the results presented in Table 6. It is interesting to observe that the r.m.s. radius, which in principle is an observable, is much less sensitive to the choice of a than de wave-function probabilities. Furthermore, the large to very large r.m.s. radii in the various situations are hardly surprising, in view of the small binding energies and the resulting very long tails of the D 0 D * 0 wave-function components (see Fig. 2

above).
Using Eq. (18), we now also evaluate the S-wave scattering length a S = − lim E→0 [k(E) cot δ 0 (E)] −1 . In case X and for a = 2.0 GeV −1 we thus find a S = 11.55 fm, which is large yet of the expected order of magnitude for a BE of 0.16 MeV. For even smaller BEs, the scattering length will further increase, roughly like ∝ 1/ √ BE. Let us here quote from Ref. [17]: "Low-energy universality implies that as the scattering length a increases, the probabilities for states other than D 0D * 0 orD 0 D * 0 decrease as 1/a . . . " Indeed, we verify from our Table 5 that -very roughlythe cc probability decreases as ∝ √ BE, and so like ∝ 1/a S .

Stability of results and nature of poles
In this section we are going to study the stability of our results, as well as the nature of the found solutions. So let us vary the two usually fixed parameters, viz. ω and m c , in such a way that the bare 1 3 P 1 mass remains unaltered at 3599 MeV, whereas that of the 2 3 P 1 changes as shown in Table. 7. Thus, in case I E 1 is lowered by 25 MeV, while in case II it rises by 100 MeV. The trajectories for these two new situations are plotted in Fig. 3. For I we observe  Table 7, and the others the standard case of Fig. 1; the solid (dashed) lines stand for normal (belowthreshold) resonances. All trajectories lie on the second Riemann sheet. The pole positions for the g values in Table 7 are marked by * . that, just as in the standard case depicted in Fig. 1, two poles are found relatively close to the real axis, of a dynamical and a confinement origin, respectively. However, now it is the 2 3 P 1 confinement pole that moves steadily downwards and settles on the real axis below threshold, whereas the dynamical pole moves to higher energies and eventually approaches the real axis. So the poles interchange their roles when going from the standard case to case I. Nevertheless, the values of g needed to get a bound state at 3871.68 MeV are not very different in the two cases, viz. 1.172 vs. 1.034. Such a behavior was already observed almost a decade ago, namely for D * s0 (2317) [6] charmed-strange meson. In a first, two-channel model calculation [27] the D * s0 (2317) showed up as a dynamical resonance, settling below the S-wave DK threshold, whereas the 1 3 P 0 cs state turned out to move to higher energies, with a large width, similarly to the standard X(3872) case in Figs. 1 and 3 above. However, in a more complete, multichannel approach [31] the situations got reversed, just as in the present case I. Also in our previous study [14] of the X(3872), with nine coupled channels, we reproduced the meson as a confinement pole. What appears to happen in the present case I is that shifting the bare 2 3 P 1 state to somewhat lower energies is just enough to deflect the confinement pole to the left and not the right when approaching the continuum pole. Clearly, there will be an intermediate situation for which the left/right deflection will hinge upon only marginal changes in the parameters, but resulting in two completely different trajectories. Therefore, identifying one pole as dynamical and the other as linked to a confinement state is entirely arbitrary, the whole system being dynamical because of unquenching. At the end of the day, the only thing that really counts is where the poles end up for the final parameters. The trajectories themselves are not observable and only serve as an illustration how a coupled-channel model as the one employed here mimics the physical situation. Suffice it to say that the lower pole, representing the X(3872), is quite stable with respect to variations in the parameters, owing to its proximity to the only and most relevant OZI-allowed decay channel. The higher pole, on the other hand, should not be taken at face value, since a more realistic calculation should include other important decay channels, such as D * D * , with threshold just above 4 GeV.
Concerning the other scenario with changed parameters, labeled II in Table 7 and depicted in the lower graph of Fig. 3, we see that the trajectories do not change qualitatively when going from the standard case to II. There is a displacement of the right-hand branch, about 100 MeV to the right on average, in accordance with the same shift of the bare 2 3 P 1 state. But the change in the lower, dynamical branch, is much less significant, though the value of g needed to produce a bound state at 3871-68 MeV now increases to 1.572 (see Table 7). We also notice from Fig. 3 that the two pole-trajectory branches hardly move towards one another, signaling less attraction between the poles due to a larger initial separation.
Inspecting again Table 7 as for the cc probability in cases I and II compared to the standard situation, we observe an increased value for case I and a decreased one for II. This is logical, since in case I the bare 2 3 P 1 state lies closer to the X(3872), whereas in case II it lies farther away. Nevertheless, the difference in cc probability between I and II is only about 3%, i.e., less than the Fig. 4. Normalized two-component radial wave function R(r), for cases II and "standard", corresponding to parameters in Table 7. Bold curves refer to case II, normal curve to Rc for standard case. Note: R f is indistinguishable within graphical accuracy for the two cases. R c R f II variation with a in the standard case. These comparisons lend further support to the stability of our results. Finally, in Fig. 4 we compare the wave function for case II with the standard one. We see there is no visible change in the R f component. As for R c , the first maximum gets somewhat reduced, but the secondary, negative bump even becomes a bit larger, owing to an inward shift of the node, lying now at about 3 GeV −1 . Yet, also in case II the cc component is still very significant, despite the large separation of more than 200 MeV between the X(3872) bound state and the bare 2 3 P 1 state.
From the latter and all previous results we may safely conclude that the cc component of the X(3872) wave function remains nonnegligible in a variety of scenarios, being even of comparable size as the D 0 D * 0 component in the inner region, save at very short distances.

Summary and conclusions
In the present paper, we have employed a simple and solvable Schrödinger model to study the wave function of the X(3872) meson, by treating it as a coupled cc-D 0 D * 0 system with J P C = 1 ++ quantum numbers. Transitions between the two channels are described with the 3 P 0 mechanism, through string breaking at a sharp distance a. The exact solutions to the equations allow us to easily study the trajectories of S-matrix poles as a function of the decay coupling constant g. Thus, a dynamical pole is found, becoming a bound state just below the D 0 D * 0 threshold, for different string-breaking distances a, and an appropriate coupling g. On the other hand, the pole arising from the bare 2 3 P 1 confinement state moves to higher energies and acquires a large imaginary part. However, the latter pole may not be very relevant physically, because of neglected additional meson-meson channels which will become important in that energy region.
As for the X(3872) radial wave function, the cc component R c turns out to be of significant size as compared to the D 0 D * 0 component R f , especially for intermediate r values. Moreover, even for other trial BEs, the global shape of R c and its relative magnitude vis-à-vis R f in the central region is remarkably stable. But the corresponding cc probability is relatively low, due to the very long tail of the D 0 D * 0 wave function at small binding. These results are along the lines of the analysis based on general arguments presented in Ref. [32]. Quantitatively, for the average [6] X(3872) binding of 0.16 MeV, a cc probability of 7.5-11.2% is found, for a in the range 0.4-0.6 fm, which is compatible with other recent approaches [29,30]. The corresponding r.m.s. radius turns out to be quite stable at about 7.8 fm, for the latter range of a values, while the S-wave scattering length of 11.6 fm, for a ≈ 0.4 fm, is in agreement with expectations for a BE of 0.16 MeV.
Finally, we have studied the nature of the S-matrix pole giving rise to X(3872), by varying some of the otherwise fixed parameters. Thus, a drastic modification of pole trajectories is observed, for relatively small parameter variations, making the X(3872) pole transform from a dynamical pole into one directly connected to the 2 3 P 1 bare confinement state. However, the corresponding changes in the cc probability and r.m.s. radius, as well as the coupling g needed to reproduce X(3872), are quite modest.
In conclusion, we should revisit the claims about X(3872) made in Ref. [16], quoted in the Introduction above, namely about the inevitability of X(3872) being a charm-meson molecule or virtual state, independently of the mechanism generating the state. Now, it is true that our analysis has confirmed some of the quantitative predictions in Ref. [16], viz. concerning the vanishing probability of wave-function components other than D 0 D * 0 as the BE approaches zero, and the related behavior of the D 0 D * 0 scattering length. However, we have also shown that the cc component is certainly not negligible and quite stable, in a variety of scenarios. Especially in electromagnetic processes, the prominence of this component for relatively small as well as intermediate r values will no doubt result in a significant contribution to the amplitudes. Moreover, as already mentioned above, the very unquenching of a 2 3 P 1 cc state will not only introduce meson-meson components into the wave function, but also a contribution of the 1 3 P 1 cc state, which can change predictions of electromagnetic transition rates very considerably [28]. We intend to study such processes for X(3872) in future work, on the basis of a model as the one used in the present paper, by employing the formalism developed and successfully applied in Ref. [33]. However, in a detailed and predictive calculation of electromagnetic X(3872) decays, the inclusion of the charged D ± D * ∓ channel will be indispensable [34].
For all these reasons, we do not consider X(3872) a charm-meson molecule, but rather a very strongly unitarized charmonium state. As a matter of fact, we do not believe any non-exotic mesonic resonance -whatever its origin -qualifies as a true meson-meson molecule, simply because such a state will inexorably mix with the nearest qq states having the same quantum numbers. Indeed, we have demonstrated above that, even with a bare cc state 200 MeV higher in mass, the resulting cc component in the wave function is still appreciable. So let us conclude this discussion by quoting and fully endorsing the following statement from Ref. [16]: "Any model of the X(3872) that does not take into account its strong coupling to charm meson scattering states should not be taken seriously." One of us (G.R.) is indebted to R. M. Woloshyn for the invitation to a most stimulating miniworkshop at TRIUMF last year, where the topic leading to the present paper was debated. Thanks are also due to E. Braaten and K. K. Seth for very useful discussions on X(3872). This work was partially supported by the Fundação para a Ciência e a Tecnologia of the Ministério da Educação e Ciência of Portugal, under contract no. CERN/FP/123576/2011. with the 1 × 1 S-matrix simply given by Real or complex poles can then be searched for numerically, by using Newton's method to find the energies for which cot δ l f (E) = i, on the appropriate Riemann sheet.

B Special functions, numerical methods, and kinematics
The confluent hypergeometric functions Φ and Ψ introduced in Appendix A are defined in Ref. [35], Eqs. (6.1.1) and (6.5.7), respectively. Thus, the function Φ is easily programmed as a rapidly converging power series, while the definition (6.5.7) of Ψ in terms of Φ and the gamma function Γ then also allows straightforward computation, by employing Gauss's multiplication formula for Γ (−ν) (see Ref. [36], Eq. (6.1.20)) so as to map the argument −ν to lying well inside the unit circle in the complex plane, whereafter a very fast converging power-series expansion of 1/Γ (−ν) (see Ref. [36], Eq. (6.1.34)) can be applied. The integrals for wave-function normalization and computation of r.m.s. radii are carried out by simple Gauss integration, choosing increasing numbers of points on a finite interval for the cc channel, and an infinite one for D 0 D * 0 . Note that, in the former case, the wave function falls off fast enough to allow convergence for a finite cutoff, whereas in the latter a suitable logarithmic mapping is used. In both cases though, because of the wave-function cusp at r = a and in order to avoid numerical instabilities, the domain of integration is split into two pieces, with up to 16 Gauss points in the inner region and 64 in the outer one, thus resulting in a very high precision of the results.
Although the X(3872) bound state can reasonably be considered a nonrelativistic system, we still use relativistic kinematics in the D 0 D * 0 channel, since parts of the resonance-pole trajectories involve relatively large (complex) momenta. For consistency, the same is done for all energies. The manifest unitarity of the S matrix is not