Unquenched quark-model calculation of X(3872) electromagnetic decays

A recent quark-model description of X(3872) as an unquenched $2\,{}^{3\!}P_1$ $c\bar{c}$ state is generalised by now including all relevant meson-meson configurations, in order to calculate the widths of the experimentally observed electromagnetic decays $X(3872) \to \gamma J/\psi$ and $X(3872) \to \gamma\psi(2S)$. Interestingly, the inclusion of additional two-meson channels, most importantly $D^\pm D^{\star\mp}$, leads to a sizeable increase of the $c\bar{c}$ probability in the total wave function, although the $D^0\bar{D}^{\star0}$ component remains the dominant one. As for the electromagnetic decays, unquenching strongly reduces the $\gamma\psi(2S)$ decay rate yet even more sharply enhances the $\gamma J/\psi$ rate, resulting in a decay ratio compatible with one experimental observation but in slight disagreement with two others. Nevertheless, the results show a dramatic improvement as compared to a quenched calculation with the same confinement force and parameters. Concretely, we obtain $\Gamma ( X(3872) \to \gamma \psi(2S) )=28.9$ keV and $\Gamma ( X(3872) \to \gamma J/\psi)=24.7$ keV, with branching ratio $\mathcal{R}_{\gamma\psi}=1.17$.


Introduction
Since its discovery [1] by the Belle Collaboration in 2003, the very narrow axial-vector [2] charmonium state X(3872) has become one of the favourite theoretical laboratories for meson spectroscopists, because of its remarkable closeness to the D 0D 0 (orD 0 D 0 ) and ρ 0 J/ψ thresholds, besides its seemingly too low mass for mainstream quark models. The now established [3] J P C = 1 ++ quantum numbers seem to imply that X(3872) is either the still unconfirmed [2] χ c1 (2 3 P 1 cc) meson, or an axial-vector charmoniumlike state of a different kind. For a recent review, see e.g. Ref. [4].
However, in order to understand the true nature of X(3872), one can ignore neither the presence of relatively nearby 1 ++ states in the theoretical charmonium spectrum, nor their strong coupling to the S-wave threshold D 0D 0 . In this spirit, the properties of X(3872) were recently studied in Refs. [5,6], by modelling it as an unquenched χ c2 state with additional meson-meson (MM) components, most importantly D 0 D 0 . 1 In the former paper [5], a momentum-space calculation of X(3872) was carried out, employing the Resonance-Spectrum Expansion (RSE), with all relevant two-meson channels included. This work showed that the hadronic decays of X(3872) can thus be described quite accurately, dispensing with ad hoc tetraquark or molecular approaches. On the other hand, the latter paper [6] focused on the X(3872) wave function, using instead a coordinate-space model and with only two channels, viz. cc and D 0 D 0 . The purpose was to study whether the charm-anticharm component would remain substantial, despite the very long tail of the D 0 D 0 component due to the small binding of less than 0.2 MeV [2]. Indeed, a cc probability of about 7.5% was found andeven more importantly -a corresponding wave-function component in the inner region of the same order of magnitude as that of the D 0 D 0 channel, thus ruling out a pure molecular scenario for X(3872). Similar interpretations of X(3872) were concluded in the unquenched model calculations of Refs. [7] and [8].
Besides the mentioned hadronic decays, X(3872) has also been observed [2] to decay in electromagnetic (EM) processes, namely to γJ/ψ and γψ(2S). Such decays are very sensitive to details of the X(3872) wave function, especially in its inner region, and so may discriminate among different microscopic models. Thus, the coordinate-space method for unquenched quarkonium states employed in Ref. [6] appears to be the indicated approach for such a calculation. Now, it was shown in Refs. [9,10] that, in a multichannel system with one almost unbound channel, more strongly bound channels should not be neglected beforehand for processes in which the wave function at short distances is important. This is of course all the more true for the D ± D ∓ channel in the X(3872) case, which is bound by only 8 MeV and so is expected to have a significant effect on the interior wave functions of the other components. Moreover, this channel contains charged mesons, which will contribute directly to EM transition amplitudes. Nevertheless, for completeness we also include all other OZI-allowed channels with combinations of pseudoscalar and/or vector charm-light as well as charm-strange mesons. As for the wave functions of the EM decay products J/ψ and ψ(2S), again all channels with pairs of groundstate D, D , D s , and D s mesons will be accounted for, since they all may develop non-negligible components.
In the present paper, we shall closely follow the formalism for EM decays of unquenched, "unitarised" quarkonium systems as developed in Ref. [11]. The organisation is as follows. In Sec. 2, a multichannel Schrödinger equation for confined qq channels coupled to free MM channels is written down and solved analytically. Section 3 is devoted to the computation and display of the multicomponent wave functions of the charmonium states X(3872), J/ψ, and ψ(2S), using the generic solutions derived in Sec. 2. In Sec. 4 the procedure [11] for EM decays of quarkonium states with MM components in the wave function is reviewed. Section 5 presents the results for the EM decays, in comparison with several other published model calculations. A summary and some conclusions are given in Sec. 6.

Multichannel Schrödinger model
Just as in Ref. [6], we consider X(3872) a coupled charmanticharm and MM system. The transitions between the cc and MM sectors are assumed to take place via 3 P 0 qq creation/annihilation at a sharp distance a, thus mimicking string breaking. Such a transition potential can be described in coordinate space via a spherical delta function. This choice makes the coupled-channel equations analytically solvable, provided that the confining cc potential allows solutions in terms of known functions that can be analytically continued. Thus, like in Refs. [6,5], we take a harmonic oscillator (HO) with universal frequency ω. The present generalisation beyond the model of Ref. [6] amounts to the inclusion of several MM channels instead of only one.
The coupled-channel Schrödinger equation to be solved reads ĥc whereĥ qq is the quark-antiquark Hamiltonian with a confining HO potential given bŷ with µ c the reduced quark mass m c q m c q /(m c q + m c q ), and h M M is the free MM Hamiltonian Furthermore, V cj in Eq. (1) is the transition potential modelled through a spherical delta shell with radius a where g cj is the relative coupling constant of the confined qq channel c to the j-th MM channel and λ is an overall coupling. The radial wave functions u c and v j result from the separation of the total wave function in spherical coordinates, i.e., ψ = u(r) r Y lm (θ, ϕ). We now first solve the equation analytically for r < a and r > a ignoring the delta shell, and then match both solutions and their derivatives at r = a, accounting for the delta function at this point. For the quark-antiquark components, we introduce the parameter and make the substitution u c (r) = F c (r)r 1+Lc e − 1 2 µcωr 2 . This way we find that the solutions are of the form where M and U are the Kummer and Tricommi confluent hypergeometric functions (same as the Φ and Ψ functions defined in Ref. [12] and employed in Ref. [6]). Given the properties of these functions, this guarantees that u c is a solution toĥ c qq u c = Eu c for r = a, regular at the origin, and vanishing at infinity.
For the two-meson wave function, we introduce the variable q j = ip j for each channel, with p j the corresponding relative momentum. Then we have The two-meson components v j (r) can be written as v j (r) = A j i Lj (q j r) r , r < a B j k Lj (q j r) r , r > a where i l (x) and k l (x) are the modified spherical Bessel functions of the first and third kind, respectively (cf. the functions J l and N l in Ref. [6]). The function i l is regular at the origin and divergent at infinity, whereas k l is irregular at the origin and falls off exponentially as x → ∞. This solution is valid as long as the energy of the state is below all two-meson thresholds. For convenience, we now simplify our notation, by writing M c for M (−ν c , l c + 3 2 , µ c ωa 2 ), i j for i Lj (q j a), and similarly U c and k j . In order to determine the coefficients a c , b c , A j , and B j , as well as the energy E for a given coupling λ or λ for a given E, we first use continuity of the wave function at r = a, i.e., Next we integrate the Schrödinger equation (1) from a − to a + , with infinitesimal. In doing so, first for the qq components, we obtain This set of equations is just a generalised eigensystem, with λ 2 the generalised eigenvalue and α c the generalised eigenvector, which can be written as where D is a diagonal matrix. We can use the latter equation to determine the value of λ for which there is a certain bound state with a chosen energy E. Alternatively, for a given λ, the energies of possible bound state can be found by employing Newton's method to search for zeros of the determinant Either way, the wave-function coefficients can next be calculated from the obtained α c , using Eqs. (15)(16)(17). Finally, the scale of the coefficients can be fixed by imposing the normalisation condition 3 X(3872), J/ψ, and ψ(2S) wave functions Now that we have derived the general solution for the multichannel wave functions resulting from Eq.
(1), we should focus on the specific wave functions of the axialvector (A) charmonium system X(3872) as well as the vector (V ) states J/ψ and ψ(2S), since we shall consider EM decays of the former charmonium into the latter two.
In order to account for all non-negligible meson-loop effects, we couple these systems to OZI-allowed channels containing pairs of the ground-state open-charm mesons D, D , D s and D s , being either pseudoscalar (P ) or vector (V ). Now, V states can decay into the combinations P P , P V , and V V , with odd orbital angular momentum L because of parity conservation, whereas A systems couple to channels of the P V and V V types, with even L. In Tables 1 and 2 we give the relative couplings of the different charmonia to the corresponding two-meson channels, viz. for J/ψ (or ψ(2S)) and X(3872), respectively. These Table 1. Squares of coupling coefficients for the vector charmonia J/ψ and ψ(2S). The generic notation D represents D 0 , Table 2. Squares of coupling coefficients for the axial-vector charmonium X(3872). Also see Table 1.
couplings have been calculated employing the formalism of Ref. [13], based on overlaps in an HO basis. For economy, each listed channel in these tables really represents three channels, with e.g. DD standing for D 0 D 0 , D ± D ∓ , and D + s D − s , all with the same coupling. Also note that in the V case of J/ψ and ψ(2S) two cc channels must be included, viz. 3 S 1 and 3 D 1 , giving rise to two sets of couplings.
Before we can compute the different wave functions, we must fix the model parameters, viz. ω (HO frequency), m c (charm-quark mass), λ (overall coupling constant), and a (string-breaking distance). Now, the former two parameters are not really free, as they have been kept fixed at the values ω = 190 MeV and m c = 1562 MeV, determined in Ref. [14], in all subsequent work. Then, λ and a should be adjusted to the masses of J/ψ, ψ(2S), and X(3872), which can be done reasonably well, in spite of having only two parameters to fit three observables. Nevertheless, we believe that in the present calculation it is most important to have as accurate as possible wave functions, so that we somewhat relax the usual condition of only one λ for all described systems. This way we are able to precisely reproduce the experimental J/ψ, ψ(2S), and X(3872) masses, with the values λ ψ = 2.527, λ X = 2.176, and a = 1.95 GeV −1 , the coupling λ ψ being of course the same for J/ψ and ψ(2S), as we have included exactly the same MM channels for these two V states. Again, ideally there should be only one λ. However, one must realise that the completeness property of the couplings g i as computed in the scheme of Ref. [13], which implies i g 2 i = 1 for a system with any quantum numbers, is only satisfied if all decay channels are included. Well, in the latter HO-based formalism, the number of allowed channels is finite but still huge and so too large to totally account for in practical calculations. Therefore, when coupling only to the most important channels, generally those with the lowest thresholds, somewhat different values of λ for clearly distinct systems are perfectly acceptable. And indeed, λ ψ and λ X differ by only about 15%. Moreover, λ X = 2.176 is of the same order of magnitude as λ ≈ 3 obtained in the X(3872) study of Ref. [5], in which the related yet quantitatively different momentum-space RSE formalism was employed, for a similar value of the decay radius, viz. a = 2 GeV −1 . Now we are in a position to compute the three needed radial wave functions, using Eq. (19) and Eqs. (15)(16)(17). In Fig. 1 the J/ψ and ψ(2S) wave functions are displayed, and in Fig. 2 that of X(3872). The first thing we observe in all three plots is a kink in the wave-function components at r = a, which is a direct consequence of our choosing a singular transition potential, mimicking string breaking at that precise distance. Concerning the V charmonia of Fig. 1, there is a dominant l = 0 cc wave function, which does not vanish at the origin, but also considerable l = 2 cc and L = 1 MM components. In the case of X(3872), which has a seemingly dominant l = 1 cc wave function, the L = 0 D ± D ∓ and D 0 D 0 components are also very sizable, especially the latter. As a matter of fact, the D 0 D 0 channel turns out to be the most important one in terms of probability, due to its very long tail, resulting from the small binding with respect to the X(3872) mass [5]. This and all other wave-function probabilities are given in Ta-  Table 1 and text). By convention and for clarity purposes, the MM curves take negative values.  Fig. 1 and Table 2. Table 3. Percentage probabilities of J/ψ and ψ(2S) wavefunction components, with the cc and MM numbers including l = 0, 2 and L = 1, 3 contributions, respectively. Also, D stands generically for D 0 , D ± , or D ± s , etc., as in Fig. 1 bles 3 and 4, i.e., for the V charmonia and X(3872), respectively. What may look surprising in Table 3 is that J/ψ has clearly larger MM components in its wave function than ψ(2S). But this can be understood by observing that the decay radius a is relatively close to the node in the ψ(2S) cc wave function, which reduces the influence of the MM channels. As for X(3872), we see in Table 4 that the D 0 D 0 component is by far the largest, just like in Table 4. Percentage probabilities of X(3872) wave-function components, with the MM figures including L = 0, 2 contributions. Also, the D D value accounts for the D 0 D 0 , the two-channel model study of Ref. [6]. Nevertheless, we observe a significant decrease in the latter MM channel's probability, and a large increase in the cc probability, viz. from about 7.5% to almost 27%, for a similar decay radius a. At first sight, this seems very surprising, as in the present work we include several additional MM channels. However, it must be realised that all such channels contribute to shift the X(3872) bare energy level of 3929 MeV down to the physical value of 3871.69 MeV, which reduces the relative importance of the D 0 D 0 channel. Since it is precisely the latter wave-function component that has a very large extension in coordinate space, the reduction of its coupling will reduce its probability roughly proportionally, mostly to the benefit of the cc component. Another way to look at the different wave-function components is by focusing on how the quenched solutions are modified by the coupling to the MM channels. For that purpose, we compute the overlaps of the unquenched cc wave functions, which are the solid curves in Figs. 1 (top, bottom) and 2, with pure HO solutions of different radial quantum number n, the results being given in Tables 5,  6, and 7, respectively. Note that these numbers do not Table 5.  concern the total wave-function overlaps, which can easily be obtained through multiplication by the cc probabilities in Tables 3 and 4. It is interesting to see that the overlaps with the lowest four radial HO states, including of course also all the l states that couple, are almost sufficient to reconstruct the cc wave functions, namely to 97% for J/ψ, and to even more than 99% for ψ(2S) and X(3872). This might be useful for quark models with confinement mechanisms different from ours, in which sometimes HO wave functions are used to compute certain observables. In the specific case of X(3872), we get an 8.82% overlap with the 1 3 P 1 bare HO state. This is interesting, as in Ref. [18] such a component, despite being smaller in size than our result, was found to have a large influence on the EM decay widths of X(3872).
To conclude this section about wave functions, we compute the r.m.s. radii of J/ψ, ψ(2S), and X(3872), obtaining 0.456 fm, 0.930 fm, and 6.57 fm, respectively. Notice that the latter number is somewhat smaller than the value found in Ref. [6], which is logical in view of the here reduced influence of the long D 0 D 0 tail.

Electromagnetic transitions
In this section we review the formalism [11] for EM transitions of quarkonium systems coupled to MM channels. In order to calculate the EM decay rate of a multicomponent meson state, we couple the EM field to our coupledchannel strong-interaction HamiltonianĤ qq−M M , obtaining a Hamiltonian of the typê whereĤ em is the free EM part, in Gaussian units readinĝ andĤ int describes the interaction between the hadrons and the EM field. This interaction Hamiltonian can be naturally obtained fromĤ qq−M M via a minimal-coupling prescription. As we know, the hadronic coupled-channel HamiltonianĤ qq−M M has diagonal elements of the form The Hamiltonian for hadrons interacting with radiation is then obtained through the minimal coupling For the EM radiation field, we use the gauge conditions Using Eqs. (25) and (26) and allowing for a possible anomalous magnetic moment µ i , we get (27) Since we are considering the EM interaction perturbatively, the term quadratic in A can and will be neglected. Now, for quarks and antiquarks, the magnetic moment is given by Then, the meson magnetic moment is obtained by assuming a pure qq state, resulting in the expression Note that, apart from accounting for the mesons' magnetic moments, the present calculation neglects their internal structure. As a consequence, to describe the EM decays of unitarised mesons, we shall only consider processes of the type (QQ) * → QQ + γ (30) and while neglecting the ones that change the internal structure of individual mesons, viz.
As the wave function of a unitarised meson has the form the total matrix element for an EM transition is given by The quantised EM vector potential A can be written in Gaussian units as [11] A(r, t) = √ 4π c λlm dk 2π where the index λ ∈ {e, m} indicates whether the component is an electric multipole or a magnetic multipole. For instance, the components with λ = e and l = 2 correspond to electric quadrupole (E2) radiation, while those with λ = m and l = 1 stand for magnetic dipole (M1) radiation. Furthermore, a λlm (k) and a † λlm (k) are photon annihilation and creation operators obeying the commutation relations The vector fields f (λ) kJM (r) are given by [11] f (m) and They have the properties [11] ∇ × f (e) Now, we have an initial state |Ψ i = |nJM ⊗ |0 , where |0 denotes the photon vacuum, and a final state Substituting next Eqs. (35), (37), and (38) into the matrix element we finally obtain, after a laborious yet straightforward calculation, the electric and magnetic decay matrix elements [11]. The results are given in Appendix A.

Results and comparison of EM decay rates
Now we can present the results for our model calculation of the X(3872) EM decays. First, though, we show in Table 8 the up-to-date experimental status of such decays, which is clearly poor. First of all, only lower bounds are reported for the γJ/ψ and γψ(2S) rates. On the other hand, for the very small total X(3872) width, merely an upper bound of 1.2 MeV is listed [2]. This is understandable, as small enough bin sizes to pin down the width more precisely are not possible with present-day statistics. But as a consequence, the absolute magnitudes of the two observed EM decays are largely unknown. Only their ratio has been determined by two experiments, though still with large errors (see Table 8). Coming now to our model predictions, we first observe that a process of the type 3 P 1 → γ 3 S 1 / 3 D 1 is dominated by an electric dipole (E1) transition, besides a smaller magnetic quadrupole (M2) contribution. Using the expressions in Appendix A, we then obtain the results presented in Table 9. As expected, the M2 widths are much smaller than the E1 ones, since higher multipoles are suppressed by powers of photon momentum divided by (charm) quark mass [30]. Such a behaviour is roughly confirmed by our numbers, as the photon in the process with J/ψ in the final state has about four times as much momentum as in the decay to ψ(2S) [2]. Clearly, experimental statistics is insufficient so far to do an angular analysis needed [31] for disentangling the E1 and M2 contributions in the X(3872) EM data. Therefore, in order to compare our prediction for the EM branchingrate ratio defined in Eq. (45) with the experimental values given in Table 8, we simply sum the E1 and M2 contributions, obtaining R γψ = 1.17. This number is compatible with the Belle [16] upper bound of 2.1, but somewhat too small as compared to the BaBar [15] and LHCb [17] measurements. However, the experimental values are only marginally in agreement with one another and the uncertainties are still very large. If we take our absolute EM-width predictions at face value, in particular the result 28.9 keV for Γ (X(3872) → γψ(2S)), the experimental lower bound B γψ(2S) > 0.030 [15] implies Γ (X(3872)) < 1 MeV, slightly lower than the PDG [2] upper limit of 1.2 MeV. Nevertheless, the X(3872) total width is probably even smaller than that, considering the S-matrix pole trajectories of the X(3872) resonance near the D 0 D 0 threshold in the multichannel RSE calculation of Ref. [5]. Next, in Table 10 we compare the present EM results to those of a number of other model calculations. Clearly, our values are somewhere in the middle of the ballpark of often disparate numbers. Generally, one may conclude that the experimental EM rate ratio seems to favour models based on a 2 3 P 1 cc assignment for X(3872), with or without other components.
Finally, we have a look at the importance of unquenching in our model. To that end, we compute the EM predictions from bare HO cc wave functions only, with unchanged parameters except for the overall couplings to MM channels, which are set to zero. Note that this results in bare J/ψ, ψ(2S), and X(3872) masses that are roughly 300, 100, and 100 MeV larger than the physical ones, respectively. Having this proviso in mind, we obtain for the total EM widths Γ quenched (X(3872) → γJ/ψ) = 0.61 keV and Γ quenched (X(3872) → γψ(2S)) = 159 keV. Comparing to the unquenched total (E1+M2) widths in Table 9, we see a rather dramatic importance of unquenching.
Concluding our study of multichannel effects, we calculate the separate contributions from the cc and the MM channels. Note that this is done in the full model and so using the wave functions plotted in Figs. 1 and 2. The total exclusive cc and MM results are then obtained by setting the wave functions of the MM and cc components to zero by hand, respectively. Thus, we get Γ cc γJ/ψ = 15.3 keV, Γ M M γJ/ψ = 1.12 keV, Γ cc γψ(2S) = 28.0 keV, and Γ M M γψ(2S) = 0.01 keV. These results reveal constructive interference effects, especially in the J/ψ case, as one cannot just sum the cc and MM numbers to obtain the widths in Table 9.

Summary and conclusions
This paper is the third part of a triptych aimed at understanding the microscopic dynamics of the enigmatic char-monium state X(3872). The first article [5] was devoted to its hadronic decays, including the OZI-forbidden ones, arriving at a good description of the available data. The second paper [6] focused on the importance of X(3872)'s cc component, concluding that a purely molecular assignment is ruled out. In our present work, we have generalised the r-space method employed in the latter paper so as to include all virtual open-charm MM channels that may contribute significantly to the X(3872) wave function, as well as to those of the vector charmonia J/ψ and ψ(2S). These wave functions have then been used to compute the widths of the experimentally observed EM decay processes X(3872) → γJ/ψ and X(3872) → γψ(2S), employing the formalism developed and applied to charmonium and bottomonium in Ref. [11].
In the first place, and concerning the thus obtained wave functions, it is quite remarkable that the inclusion of several additional MM channels to describe X(3872), most notably the D ± D ∓ channel bound by only about 8 MeV, gives rise to an increase of the cc probability from 7.48% [6] to 26.76%. Accordingly, the cc wave-function component becomes clearly the dominant one, except at very small and very large distances. However surprising at first sight, this effect can be explained by the reduced influence of the narrowly bound D 0 D 0 channel, which has an extremely long tail and so takes up most of the total probability. This is also confirmed by the reduction of the X(3872) r.m.s. radius from 7.82 fm in Ref. [6] to 6.57 fm here. As for the J/ψ and ψ(2S) wave functions, we find significant 3 D 1 cc and P -wave MM components, the former most prominently for ψ(2S) and the latter especially in the J/ψ case. The closeness of the ψ(2S) wave-function node to the decay radius r offers an explanation for the relatively reduced importance of MM channels for the first radially excited vector charmonium state.
Coming now to the EM transitions of X(3872), we obtain total widths of 28.9 keV and 24.7 keV for the decays to γψ(2S) and γJ/ψ, respectively, with rate ratio R γψ = 1.17. While there are no data for the absolute magnitudes of the EM widths, the latter ratio can be compared to three experimental values. Thus, our prediction of 1.17 is fully compatible with the upper bound of 2.1 observed by the Belle Collaboration [16], but a little bit too low for the numbers 3.4 ± 1.4 and 2.46 ± 0.64 ± 0.29 reported by BaBar [15] and LHCb [17], respectively. Although the experimental error bars are still quite large and the three observations so far are nonetheless only in marginal agreement with one another, additional mechanisms -not considered here -might remove the slight discrepancy. One possibility is the inclusion of photonic decays of individual mesons in the MM channels, that is, observed [2] processes of the type D + → D + γ, D 0 → D 0 γ, and D + s → D + s γ. Another source may be relativity, as relativistic effects are capable of shifting the nodes in the X(3872) and ψ(2S) cc wave-function components [11], thus affecting the overlap integrals. While these issues are certainly worthwhile to be studied in future work, the most important contribution to an even better understanding of X(3872) will be improved measurements, with higher statistics and smaller bin sizes, in order to pin down the absolute magnitudes of the EM decay widths, as well as the total width.

A EM transition matrix elements
The multipolar electric transitions matrix elements are given by [11] with the radial integrals The multipolar magnetic transitions matrix elements are [11] M m if = i Accounting for the recoil of the final-state meson, the photon momentum k is given by The decay width is then given by the Fermi Golden Rule where ρ f is the density of final states ρ f = 1/2π c.