O ct 2 02 0 Analysis of the γγ → DD̄ reaction and the DD̄ bound state

En Wang, 2, ∗ Hong-Shen Li, Wei-Hong Liang, 3, † and Eulogio Oset 4, ‡ School of Physics and Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China Department of Physics, Guangxi Normal University, Guilin 541004, China Guangxi Key Laboratory of Nuclear Physics and Technology, Guangxi Normal University, Guilin 541004, China Departamento de F́ısica Teórica and IFIC, Centro Mixto Universidad de Valencia CSIC, Institutos de Investigación de Paterna, Aptdo. 22085, 46071 Valencia, Spain


I. INTRODUCTION
The χ c0 (2P ) state was introduced in the PDG [1] (χ c0 (3860)) with J P C = 0 ++ , M = 3862 +26+40 −32−13 MeV, and Γ = 201 +154+88 −67 −82 MeV] 1 based on the single experimental measurement on the e + e − → J/ψDD reaction reported by the Belle Collaboration [2], by looking into the DD invariant mass distribution close to threshold. These experimental data have only four points below 3900 MeV, with very large errors [2], and it was shown in Ref. [3] that this information was not sufficient to draw any conclusion about the existence of this state. Three important facts were stressed in Ref. [3]: 1) The data divided by phase space did not show any peak that would justify a claim of a state at 3860 MeV. 2) A fit of the data with a bound state of DD, which has been found in Refs. [4][5][6], was possible, but again the uncertainties were too large to make any conclusive claims. 3) A fit of the data close to threshold using a Breit-Wigner, as discussed in Ref. [3], should be avoided. This last point has been often recalled concerning fits to data [7,8].
The question remains whether there are other data which can provide good information on the possible DD bound state. One attempt was done in Ref. [9] using early data of the e + e − → J/ψDD reaction from the Belle Collaboration [10]. Although a DD bound state was found consistent with the data, the quality of these data did not allow one to be too strong on the claim of this bound state. Several reactions have been suggested, measuring DD mass distributions close to threshold, which can help, with good statistics, to bring an answer to this question. In Ref. [11] three methods were devised to find an answer to this problem; the first one is the radiative decay of the ψ(3770), ψ(3770) → γX(3700) → γηη. The second one proposes the analogous reaction ψ(4040) → γX(3700) → γηη, and the third reaction is the e + e − → J/ψX(3700) → J/ψηη. In Ref. [12] the B 0 decay to D 0D0 K 0 reaction was suggested. The B + → D 0D0 K + reaction has been measured by the BaBar Collaboration [13] and is well reproduced in Ref. [12], but the unmeasured B 0 → D 0D0 K 0 reaction was found to be more useful because it does not have the tree level contribution for D 0D0 production and is proportional to the D + D − → D 0D0 transition amplitude which contains the bound state. In Ref. [14] the ψ(3770) → γD 0D0 decay was retaken, separating the D + D − production from the D 0D0 one and showing that the latter has a much bigger potential to provide valuable information concerning the existence of the DD bound state.
Awaiting future results from some of the suggested reactions, there are interesting data that we wish to investigate here concerning that point, and these are the γγ → DD data measured by the Belle [15] and BaBar Collaborations [16]. In Ref. [15], the Belle Collaboration has reported the charmonium state X(3930) in the reaction of γγ → DD, with mass 3929 ± 5 ± 2 MeV and width 29 ± 10 ± 2 MeV, which are consistent with expectations for the χ c2 (2P ) charmonium state. Later the BaBar collaboration has also performed the γγ production of the DD system, and the DD invariant mass distribution shows clear evidence for the X(3930) state, its mass and width determined to be M = 3926.7 ± 2.7 ± 1.1 MeV, and Γ = 21.3 ± 6.8 ± 3.6 MeV [16].
On the other hand the Belle and BaBar data of Refs. [15,16] were also used in Ref. [17], making fits with Breit-Wigner structures, to suggest that there could be an indication of a χ c0 (2P ) state around 3840 MeV and a width about 200 MeV, with the warning that "More refined analysis of the data with higher statistics is definitely necessary to confirm our assertion". Actually we will show that the data divided by phase space does not show any peak around 3840 MeV, thus weakening the guess of Ref. [17], and then provide an alternative explanation of the combined data of Belle and BaBar [15,16] based on the explicit consideration of the DD final state interaction, which can shed some light on the possible DD bound state.

II. FORMALISM
A. DD interaction in I = 0 In Ref. [4], the s-wave meson-meson scattering in the charm sector was studied and a prediction for a DD bound state with isospin I = 0 was made. In Ref [9], it was found that the state with mass M DD = 3730 MeV, and width Γ DD = 30 MeV was compatible with the data of the process e + e − → J/ψDD reported by the Belle Collaboration [10]. In Ref. [11], where three methods to detect this state were suggested, a state with M DD = 3720 MeV and Γ DD = 36 MeV was found including decays to all possible pairs of light pseudoscalar.
In this paper, we use only one channel, apart from the three channels D + D − , D 0D0 , and D sDs , which is ηη to account for the width of the DD bound state, as used in Refs. [3,11,12]. The transition potentials V i,j (i, j = D + D − , D 0D0 , and D sDs ) are tabulated in Table 9 of Appendix A of Ref. [4], and we introduce the potentials of ηη → D + D − and ηη → D 0D0 with a dimensionless strength a ηη to give the width of the DD bound state. The transition potentials of ηη to ηη and D sDs are not relevant and are taken as zero. As done in Ref. [3], we will multiply the potentials V D + D − ,DsDs and V D 0D0 ,DsDs by a factor f DsDs to stress more the cusp effect.
Then the amplitude t i,j for the i channel to j channel can be obtained from the Bethe-Salpeter equation, where the matrix G is diagonal with each of its elements given by the loop function for the two particles, and we take the expression of the dimensional regularization as shown in Eq. (31) of Ref. [4], where the µ = 1500 MeV, and the subtraction constant α will be taken as a free parameter.

B. Model for the reaction γγ → DD
In this section, we present the model for the reaction, where p, k, p ′ , and k ′ are the four-momenta of the two incoming photons, D + , and D − , respectively, and ǫ 1,2 are the polarizations of the two incoming photons. We can get the mechanism for this process inspired by the work of Ref. [18], where the reactions γγ → π + π − , π 0 π 0 were studied. In Ref. [18], the whole range of ππ invariant mass from 280 MeV (2m π ) till 1400 MeV was studied. The model used was good up to about 1000 MeV with no free parameters, and for higher energies the f 2 (1270) excitation was introduced by hand. In the present case, we only need the model for about a range of 144 MeV, from the DD threshold to 3880 MeV. The model for the process γγ → DD combines the Born terms: the contact term, the D meson exchange in the t and u channels, as shown in Fig. 1. Contrary to the γγ → π + π − , the D exchange terms are now here much smaller than those of π exchange of Ref. [18], because we have the denominator in the D meson propagator at threshold, where q = (q 0 , q) is the four-momentum of the exchanged D meson, and we have q 0 = p 0 − p ′0 = 0 and | q| = | p| = p 0 ≈ m D at the DD threshold. So we have, which is much smaller than 1/(−2m 2 π ) in absolute value. These terms have also energy dependence, because we have the vertex with the term ǫ · (p ′ − q), which in the Coulomb gauge ǫ 0 = 0 and ǫ · p = 0 for the photon, which we use to evaluate, is given by, In the limited range of the DD invariant masses that we consider p ′ is small and one can easily see that the contribution of the D-exchange terms are smaller than 3% of the contact term of Fig. 1(a), 2e 2 ǫ 1 · ǫ 2 . Hence we neglect these exchange terms and take the amplitude as, Thus, we will neglect the contributions of Figs. 1(b) and (c) in this work.
In addition, there are also other possible exchanges of D * resonances with anomalous terms but again, the denominator of the propagators are large and the terms are small close to the threshold.
We have the differential cross section for the reaction γγ → DD, where we average the polarization vectors of the transverse photons, with no angular dependence. Thus, we have the cross section, where s = (p + k) 2 , p and p ′ are the three-momenta of the incoming photon and the D + in the center-mass frame, respectively, In our process, we have to take into account the final state interaction of the mesons D andD. We will differentiate between D + D − or D 0D0 , since in the experiments there is information about both.

C. Final state interaction
So far we have evaluated the amplitude and cross section of the γγ → D + D − at the tree level without considering the final state interaction of D + D − . We address here this problem. In addition, the γγ → D 0D0 is null at this level and it can only proceed via rescattering of D + D − → D 0D0 . This makes this reaction more favorable to learn about a possible DD bound state since the γγ → D 0D0 amplitude is then proportional to the D + D − → D 0D0 amplitude which contains information on this possible state. The final state interaction proceeds as depicted in Fig. 2(b). The amplitude in Eq. (6) is now replaced by, where, and hence, Equations (12) and (13) can be rewritten as, The cross section is now given by Eq. (9) multiplying it by |t D + D − | 2 or |t D 0D0 | 2 for D + D − or D 0D0 production, respectively. The interpretation of the data in Refs. [15,16] requires a prior discussion. The first surprise is that in both experiments there are more events of D 0D0 production than for D + D − production. This is surprising since the strengths of t D + D − →D + D − and t D + D − →D 0D0 are the same (see Eq. (16)), but in the case of D + D − production we have the additional tree level mechanism (see Eq. (12)). The answer to this question has to be seen in Table II of Ref. [15] where the D decay modes used in the detection are shown (the same detection method is used in Ref. [16]). For D 0D0 production, four decay modes are considered: 1) D 0 → K − π + ,D 0 → K + π − ; 2) D 0 → K − π + ,D 0 → K + π − π 0 ; 3) D 0 → K − π + ,D 0 → K + π − π − π + ; 4) D 0 → K − π + π + π − ,D 0 → K + π − π 0 . However, for the D + D − production only the D + → K − π + π + , D − → K + π − π − decay mode is considered. It is thus not surprising that more D 0D0 production events than D + D − ones are observed. Inspection of the data in Fig. 5 of Ref. [15] shows that the strength of the D + D − production around 3850 MeV is about 1/3 of that of the D 0D0 production. We shall take this into account when comparing with the data. To increase the statistics, the sum of the two production modes is shown in Fig. 10 of Ref. [15], and we shall compare with those data taking into account the experimental weights for the D + D − and D 0D0 production. On the other hand, the data of Ref. [16] for D 0D0 production have a good statistics to compare directly with them. In view of that, in order to compare with the BaBar [15] and Belle [16] data, we shall use Eq. (9) multiplied by |t Belle | 2 and |t BaBar | 2 , where, with a factor B adjusted to get σ D + D − about 1/3 of σ D 0D0 around 3850 MeV. The normalization factors C Belle and C BaBar are introduced to compare with the number of events in Ref. [15] and Ref. [16] instead of cross sections.  Fig. 2(a) of Ref. [15], and the BaBar data for γγ → DD are taken from Fig. 10 of Ref. [16]. The units are in an arbitrary normalization.

III. RESULTS AND DISCUSSIONS
In this section, we will show our results. Firstly, we divide the DD invariant mass distributions of Belle and BaBar by the phase space factor | p ′ |/(s| p |) of Eq. (9), which, up to an arbitrary normalization, are shown in Fig. 3(a) and Fig. 3(b), respectively for the Belle and BaBar data. One can find that there are no peaks around 3860 MeV, and both distributions peak at the threshold, which implies that some possible states below the threshold may play an important role in the reaction of γγ → DD, and the similar feature was found in thepΛ invariant mass distribution of χ c0 →pKΛ [19].
As discussed above, there are five parameters: 1), a ηη the dimensionless potential of ηη → D + D − and ηη → D 0D0 ; 2), an extra factor f DsDs of the potentials V D + D − ,DsDs and V D 0D0 ,DsDs ; 3), the subtraction constant α in the loop function; 4), two normalization factors C Belle and C BaBar . We will fit these parameters to the experimental data in following. It should be noted that the amplitudes produced by our model have a limited range of validity and  Table I. The curves labeled 'our' stand for our calculations for the DD invariant mass distribution or the modulus squared of the amplitude, the error band labeled as '68%CL' reflects the uncertainties from the fit and represent the 68% confidence level, the Belle data are taken from Fig. 2(a) of Ref. [15]. should not be used much above the D sDs threshold, thus we only consider the experimental data points from the DD threshold to 3860 MeV.
In the first step, we fit to the Belle data of γγ → D 0D0 alone (Fit A). The fitted parameters are tabulated in Table I, and the mass distribution is shown in Fig. 4(a). Our results are in good agreement with the Belle data of γγ → D 0D0 . With the fitted parameters, the modulus squared of the amplitudes |t DD→DD | 2 is depicted in Fig. 4(b), where we can find that there is a peak around 3730 ∼ 3740 MeV, associated to a bound DD state.
Next we perform the fit to the BaBar data of γγ → DD alone (Fit B), and the fitted parameters are tabulated in Table I. With the fitted parameters, we show the DD mass distribution and the modulus squared of the amplitude |t DD→DD | 2 in Fig. 5(a) and Fig. 5(b). We have adjusted the relative weight B of Eq. (20) to get σ D + D − about 1/3 of σ D 0D0 around 3850 MeV in this case and also in the following fits. It is easy to see that there a peak around 3720 MeV, which can also be associated to the DD bound state.
Then we perform the fit to both the Belle and Babar data (Fit C), and the fitted parameters are tabulated in Table I. We present the DD invariant mass distributions in Fig. 6(a) and Fig. 6(b), respectively for the Belle and BaBar. With the fitted parameters, the modulus squared of the amplitude |t DD→DD | 2 is given in Fig. 6(c). Taking into account the uncertainties, our results are in reasonable agreement with the Belle and BaBar measurements, and the fit favors a narrow bound DD state, which can be seen from Fig. 6(c). As we discussed in Ref. [3], the present quality of the e + e − → J/ψDD data from the Belle Collaboration [2] did not allow one to be too strong on the claim of a DD bound state around 3720 MeV, although this DD bound state was found to be compatible with the Belle measurements. Since the DD final state of the e + e − → J/ψDD reaction is the same as the one of γγ → DD, we make a global fit to the data of γγ → D 0D0 of Belle [15], γγ → DD of BaBar [16], and e + e − → J/ψDD of Belle [2] 2 (Fit D), and the fitted parameters are tabulated in Table I. The mass distributions of γγ → DD are shown in Fig. 7(a) and Fig. 7(b) for Belle and BaBar, respectively. The DD mass distribution of e + e − → J/ψDD is shown in Fig. 7(c). With the fitted parameters, the modulus squared of the amplitudes |t DD→DD | 2 is given in Fig. 7(d). The global fit also favors a DD bound state around 3720 MeV.
show that there are still large uncertainties to be assertive about the position and width of the state

IV. CONCLUSIONS
In this work, we have investigated the reaction of γγ → DD by taking into account the s-wave DD final state interactions. Since the present quality of the e + e − → J/ψDD data from the Belle Collaboration did not allow one to be too strong on the claim of the DD bound state, and the final states of e + e − → J/ψDD and γγ → DD are the same, we perform four kinds of fits to the data of γγ → D 0D0 from the Belle Collaboration, γγ → DD from the BaBar Collaboration, and e + e − → J/ψDD from the Belle Collaboration. Considering the uncertainties from the fitted parameters, our results are consistent with the experimental data in the four fits, and the modulus squared of the amplitudes |t DD→DD | 2 show peaks around 3710 ∼ 3740, which can be associated to the DD bound state. Yet, the explicit evaluation of the errors done in each of the fit to the data, and particularly the last one including all the data, show that there are still large uncertainties to be assertive about the position and width of the state. Thus we encourage our experimental colleagues to measure both reactions with larger data samples, which can allow us to reduce the errors and be more precise concerning the DD bound state.