Exclusive Initial-State-Radiation Production of the $D \bar D$, $D \bar D^*$, and $D^* \bar D^*$, Systems

We perform a study of the exclusive production of $D \bar D$, $D \bar D^*$, and $D^* \bar D^*$ in initial-state-radiation events, from $e^+ e^-$ annihilations at a center-of-mass energy near 10.58 GeV, to search for charmonium and possible new resonances. The data sample corresponds to an integrated luminosity of 384 $fb^{-1}$ and was recorded by the BaBar experiment at the PEP-II storage rings. The $D \bar D$, $D \bar D^*$, and $D^* \bar D^*$ mass spectra show clear evidence of several $\psi$ resonances. However, there is no evidence for $Y(4260) \to D \bar D^*$ or $Y(4260)\to D^* \bar D^*$.


I. INTRODUCTION
The surprising discovery of new states decaying to J/ψπ + π − [1,2] has renewed interest in the field of charmonium spectroscopy, since the new resonances are not easy to accommodate in the quark model. In particular, the BABAR experiment discovered a new broad state, Y (4260), decaying to J/ψπ + π − in the initial-stateradiation (ISR) reaction e + e − → γ ISR Y (4260). The quantum numbers J P C = 1 −− are inferred from the single virtual-photon production mechanism. Further structures at 4.36 GeV/c 2 [3,4] and 4.66 GeV/c 2 [4] have been observed in the ψ(2S)π + π − mass distribution from the reaction e + e − → γ ISR ψ(2S)π + π − . Charmonium states at these masses would be expected [5,6] to decay predominantly to DD, D * D, or D * D * [7]. It is peculiar that the decay rate to the hidden charm final state J/ψπ + π − is much larger for the Y (4260) than for excited charmonium states [8], and that at the Y (4260) mass the cross section for e + e − → hadrons exhibits a local minimum [9]. Several theoretical interpretations for the Y (4260) have been proposed, including unconventional scenarios: quark-antiquark gluon hybrids [10], baryonium [11], tetraquarks [12], and hadronic molecules [13]. For a discussion and a list of references see, for example, Ref. [14].
This work explores ISR production of DD, D * D, and D * D * final states for evidence of charmonium states and unconventional structures. This follows an earlier BABAR measurement of the DD cross section [15]. A study by the Belle collaboration of the DD, D * D, and D * D * final states can be found in Refs. [16,17]. Recent measurements of the e + e − cross sections can be found in Ref. [18].
We also measure for the first time branching fractions of high mass charmonium states, other than Y (4260), for which little information exists [9], and compare our measurements with theoretical expectations [5,6,14].
This paper is organized as follows. In Section II we give a short description of the BABAR experiment and in Section III we describe the data selection. Section IV is devoted to the selection of the D * D final state and in Section V, we present the mass resolution, reconstruction efficiency, and measured cross sections. In Section VI we describe the D * D * cross section measurement while in Section VII we present the DD data. The description of the fit of the three channels is described in Section VIII, while Section IX is devoted to the measurements of the ratios of branching fractions. Finally, in Section X, we compute the limit on production of Y (4260) decaying to D * D and D * D * , and summarize conclusions in Section XI.

II. THE BABAR EXPERIMENT
This analysis is based on a 384 fb −1 data sample recorded at the Υ(4S) resonance and 40 MeV below the resonance by the BABAR detector at the PEP-II asymmetric-energy e + e − storage rings. The BABAR detector is described in detail elsewhere [19]. We mention here only the parts of the detector which are used in the present analysis. Charged particles are detected and their momenta measured with a combination of a cylindrical drift chamber (DCH) and a silicon vertex tracker (SVT), both operating within a 1.5 T magnetic field of a superconducting solenoid. The information from a ringimaging Cherenkov detector combined with energy-loss measurements in the SVT and DCH provide identification of charged kaon and pion candidates. Photon energies are measured with a CsI(Tl) electromagnetic calorimeter.
III. DATA SELECTION DD candidates are reconstructed in the seven final states listed in Table I. The D * 0 → D 0 π 0 and D * 0 → D 0 γ decay modes are used to form D * 0 D 0 and D * 0 D * 0 candidates. The D * + → D 0 π + and D * + → D + π 0 decay modes are used to form D * + D − and D * + D * − candidates. Table II summarizes the full decay chains used to reconstruct the D * D and D * D * candidates.
For all final states, events are retained if the number of well-measured charged tracks, having a minimum transverse momentum of 0.1 GeV/c, is exactly equal to the total number of charged daughter particles. Photons are identified as EMC clusters that do not have a spatial match with a charged track, and that have a minimum energy of 30 MeV. Neutral pion candidates are formed from pairs of photons kinematically fitted with the π 0 mass constraint. K 0 S candidates are reconstructed, with a vertex fit, in the π + π − decay mode. The tracks corresponding to the charged daughters of each D candidate are constrained to come from a common vertex. Additionally, for the D 0 → K − π + π 0 channel, the D 0 mass constraint is included in the fit, and for the D − → K 0 S π − channel, a K 0 S mass constraint is imposed. Reconstructed D candidates with a χ 2 fit probability greater than 0.1% are retained. Each DD pair is refit to a common vertex with the constraint that the pair originates from the e + e − interaction region. Only candidates with a χ 2 fit probability greater than 0.1% are retained. Background π 0 candidates from random combinations of photons and other background channels are suppressed by requiring no more than one π 0 candidate other than those attributed to the D 0 and D * decays. Similarly, we require in the event no more than one extra photon candidate, having a minimum energy of 100 MeV, apart from any photon attributed to D * or π 0 decays. For D decay modes without a π 0 daughter, the Dcandidate momentum is determined from the summed three-momenta of the decay particles and its energy is computed using the nominal D mass value [9]. For the D 0 → K − π + π 0 channel, the 4-momentum from the mass-constrained fit is used. Similarly, the D * momentum is determined from the summed three-momenta of the decay particles and its energy is computed using the nominal D * mass.
The ISR photon is preferentially emitted at small angles with respect to the beam axis, and escapes detection in the majority of ISR events. Consequently, the ISR photon is treated as a missing particle. We define the squared mass (M 2 rec ) recoiling against the DD, D * D, and D * D * systems using the four-momenta of the beam particles (p e ± ) and of the reconstructed D (p D ) and D * (p D * ): This quantity should peak near zero for both ISR events and for exclusive production of e + e − → D * D or e + e − → D * D * . In exclusive production the D * D and D * D * mass    Table I for the D * 0 D and D * 0 D * 0 states. The column headed "Veto" lists ambiguities with the indicated channels, "Removed" indicates the fraction of events removed by the veto.
N Channel First decay mode Second decay mode Veto Removed % distributions peak at the kinematic limit. Therefore we select ISR candidates by requiring DD, D * D and D * D * invariant masses below 6 GeV/c 2 and |M 2 rec | < 1 GeV 2 /c 4 .
We select D and D * candidates based on the candidate D mass and the mass difference ∆m = M D * − M D . The D and D * parameters are obtained by fitting the relevant mass spectra (see Fig. 1 for some ∆m distributions) using a polynomial for the background and a single Gaussian for the signal. Events are selected within ±2.5σ from the fitted central values, where σ is the Gaussian width. For D * + → D 0 π + , the selection criterion has been extended to ±6σ due to the presence of non-Gaussian tails.
Because of our tolerance of an extra π 0 and/or γ, an ambiguity can occur for channels involving a D * 0 which is handled as follows. Each combination is considered as a possible candidate for channels 8-12 and D 0 D 0 . Monte Carlo simulations weighted by the DD, D * D, and D * D * measured cross sections [15][16][17] and branching fractions are used to optimize the selection criteria and estimate the feedthrough of one channel to the other. A candidate is rejected if (a) it satisfies all the selection criteria for an ambiguous channel and (b) this rejection does not produce any significant loss in the channel under study and therefore can be classified as background. The list of channels rejected in case of ambiguities are indicated in the "Veto" column in Table II. The table also lists the fraction of events removed by these cuts in the |M 2 rec | < 1 GeV 2 /c 4 region.
In the case of multiple D * 0 candidates, such as D * 0 D * 0 with both D * 0 → D 0 γ, the candidate with m(D 0 γ) closest to the nominal D * 0 mass is accepted. The charged D * D and D * D * modes, also listed in Table II, do not require such a procedure because backgrounds are negligible.

IV. STUDY OF THE D * D FINAL STATE
The shaded regions indicate the ranges used to select the D * candidates.
gion. The resulting yields and fitted purities, defined as P = N signal /(N signal + N background ), for each channel are summarized in Table III.
The purity of each reconstructed D * channel is also demonstrated in Fig. 1, where the ∆m distribution is shown for D * D candidates with |M 2 rec | < 1 GeV 2 /c 4 and D * D masses below 6 GeV/c 2 . The final selection of the ISR candidates is performed applying the ∆m selection criteria described above.
The D * 0 D 0 mass spectrum is shown in Fig. 3(a), and the D * + D − mass spectrum is shown in Fig. 3(b). Both spectra show an enhancement near threshold due to the presence of the ψ(4040) resonance.
The background shape for D * 0 D 0 candidates is explored using the M 2 rec sideband region, 1.5 < M 2 rec < 3.5 GeV 2 /c 4 . The D * 0 D 0 mass spectrum for these events, normalized to the background estimated from the fit to the M 2 rec distribution, is presented as the shaded histogram in Fig. 3(a). This background has been fitted with a threshold function: where m th is the threshold D * 0 D 0 mass. The D * + D − final state is consistent with having zero background.
Total DD 853 77.6 ± 1.4 25 In order to measure efficiencies and D * D mass resolutions, ISR events are simulated at five different values of the D * D invariant masses between 4.25 and 6.25 GeV/c 2 . These events are simulated using the GEANT4 detector simulation package [20] and are processed through the same reconstruction and analysis chain as real events. The mass resolution is determined from the difference between generated and reconstructed D * D masses. The D * D mass resolutions are similar for all channels and increase with D * D mass from 5 to 10 MeV/c 2 .
The mass-dependent reconstruction efficiency for channel i, ǫ i (m D * D ), is parameterized by a second-order polynomial, and is multiplied by the decay branching fraction B i [9] to define These values are weighted by N i (m D * D ), the number of D * D candidates in channel i, to compute the average efficiency as a function of m D * D , where n is the number of decay modes. In this case we have eight D * 0 D 0 chanels (1-4 with D * 0 → D 0 γ and D * 0 → D 0 π 0 ) and two D * + D − channels (13,14). Representative values of ǫ B , computed at a mass of 4.5 GeV/c 2 , are displayed in Table III. The sample sizes for the Cabibbo-suppressed decay modes (15,16, and 17 in Table II) are very small (32 events total) and comprise 14% of the D * + D − sample. The efficiency for these decay channels has not been directly computed; instead, these modes are treated as having the mean efficiency of the Cabibbo-allowed channels 13 and 14. The ten D * D channels, after correcting for efficiency and branching fractions, have yields that are consistent within the statistical errors.
The D * D cross section is computed using where dN/dm D * D is the background-subtracted yield.
The differential luminosity is computed as [21] dL where s is the square of the e + e − center-of-mass energy, α is the fine-structure constant, x = 1 − m 2 D * D /s, m e is the electron mass, and L is the integrated luminosity of 384 fb −1 . The cross sections for D * 0 D 0 , D * + D − , and combined D * 0 D 0 and D * + D − are shown in Fig. 4. A clear ψ(4040) resonance is observed.
The systematic uncertainties on the cross sections, 10.9% for D * 0 D 0 and 9.3% for D * + D − , include uncertainties for particle identification, tracking, photon and π 0 reconstruction efficiencies, background estimates, branching fractions, and a potential inaccuracy in the simulation of extraneous tracks, photon and π 0 candidates. The uncertainty due to the ISR selection has been estimated by narrowing the M 2 rec allowed range to 0.7 GeV 2 /c 4 . All contributions are added in quadrature. Systematic uncertainties are summarized in Table IV.
The D * 0 D 0 and D * + D − cross sections have similar features and consistent yields. Integrating the cross sections from threshold to 6 GeV/c 2 , we obtain consistent with unity. In this calculation systematic errors related to the M 2 rec selection criteria and tracking efficiency have been ignored because they largely cancel in the ratio.

VI. STUDY OF THE D * D * SYSTEM
A similar analysis is carried out for D * D * channels. Figure 5 shows the ∆m distributions for D * D * candidates with |M 2 rec | < 1 GeV 2 /c 4 and D * D * masses below 6 GeV/c 2 . The peak at threshold in Fig. 5(a) is due to background from D * 0 → D 0 π 0 where one γ from the low momentum π 0 is lost.
We select the two D * candidates and reject candidates reconstructed in any of the modes listed in the "veto" column in Table II. Figure 6  the D * 0 D * 0 channel is estimated by fitting the M 2 rec distribution. The fit is performed using a 2 nd -order polynomial for the background and a signal M 2 rec lineshape obtained from Monte Carlo simulations that reflect the composition of the data. The number of ISR candidates and purities are also summarized in Table III. The D * + D * − final state has a background consistent with zero.
Due to the small D * + D * − sample size, the charged and neutral mass spectra are summed in Fig. 8. The D * D * mass spectrum shows unresolved peaks at ψ(4040) and ψ(4160) and an enhancement at the position of the ψ(4400) [9].
The background is explored using events in the M sideband regions −2.5 < M 2 rec < −1.5 GeV 2 /c 4 and 1.5 < M 2 rec < 2.5 GeV 2 /c 4 , and fitted using Eq. (2). The D * D * mass spectrum for these events, normalized from the fit to the M 2 rec distribution, is shown as the shaded histogram in Fig. 8.
The D * D * cross section is calculated using the same method used to compute the D * D cross section. The result, summed over the neutral and charged modes, is shown in Fig. 9. All systematic uncertainties which have been taken into account for the D * D * mode are listed in Table V; the overall uncertainty on the cross section is 12.4%.
The D * D * cross section distribution exhibits a threshold enhancement due to the superposition of the ψ(4040) The shaded regions indicate the ranges used to select the D * signals.

VII. THE DD MASS SPECTRUM
In the selection of the D 0 D 0 sample we also apply the method of resolving ambiguous events having an additional π 0 and/or γ. Here we veto all events that are ambiguous with channels 8-12, obtaining a rejection of 7.6 % background events in the |M 2 rec | < 1 GeV 2 /c 4 region. No such procedure is applied to the D + D − sample. The DD analysis is otherwise identical to that reported in Ref. [15]. The resulting M 2 rec distributions for D 0 D 0 and D + D − channels are shown in Fig. 10. The curves  are the results from the fits performed using a 2 nd order polynomial for the background and a M 2 rec lineshape obtained from Monte Carlo simulations that reflect the channel composition of the data. Again, the resulting event yields and purities are summarized in Table III. The combined DD mass spectrum is shown in Fig. 11. The background is explored using events in the M 2 rec sideband regions 1.5 < M 2 rec < 3.5 GeV 2 /c 4 and fitted using Eq. (2). This background, normalized from the fit to the M 2 rec distributions, is shown as the shaded histogram in Fig. 11. The features in the DD mass spectrum and the resulting DD cross section have been extensively discussed in our previous publication [15].  spectra are performed. We write the likelihood functions as where m is the D ( * ) D ( * ) mass, c i and φ i are free parameters, W i (m) are P-wave relativistic Breit-Wigner distributions [9], P (m) represents the nonresonant contribution, B(m) describes the background, ǫ B (m) is the av-erage efficiency, and f is the signal fraction fixed to the values obtained fitting the M 2 rec distributions. The parameters of the ψ(4040), ψ(4160), and ψ(4415) are fixed to the values reported in the Review of Particle Physics [9]. The parameters of the ψ(3770) are fixed to the values obtained in our previous analysis of the DD system [15]. The DD data require that we include the 3.9 GeV/c 2 structure, as suggested in Ref. [22], which we parametrize empirically as the square root of a Gaussian times a phase factor G(m)e iφ2 . The parameters of the Gaussian are fixed to the values obtained in our previous analysis of the DD system: m G(3900) = 3943 ± 17 MeV/c 2 , σ G(3900) = 52±8 MeV/c 2 [15]. The shape of the nonresonant contribution P (m) is unknown; we therefore parametrize it in a simple way as where C(m) is the phase space function for D ( * ) D ( * ) , and a and b are free parameters. Resolution effects have been ignored since the widths of the resonances are much larger than the experimental resolution.
Interference between the resonances and the nonresonant contribution P (m) is required to obtain a satisfactory description of the data. The size of the nonresonant production is determined by the fit.
The six D 0 D 0 , D + D − , D * 0 D 0 , D * + D − , D * 0 D * 0 , and D * + D * − likelihood functions are computed with different thresholds, efficiencies, purities, backgrounds, and numbers of contributing resonances appropriate for each channel. The fits, summed over the charged and neutral final states, are shown in Fig. 12; they provide a good description of all the data. In the figure, the shaded areas indicate the background estimated by fitting the M 2 rec sidebands. The second smooth solid line represents the nonresonant contribution where we plot |P (m)| 2 , therefore ignoring the interference effects. The fraction for each resonant contribution i is defined by the following expression: The fractions f i do not necessarily add up to 1 because of interference between amplitudes. The error for each fraction has been evaluated by propagating the full covariance matrix obtained by the fit. The resulting fit fractions and phases are given in Table VI.

IX. FIT FRACTIONS AND INTEGRATED RATES
The fit gives corrected yields for each charmonium resonance. Since the fits have been performed independently for the neutral and charged modes, the weighted means of the fit fractions are used. These can be used to compute the integrated rates for each resonance in the DD, D * D, and D * D * decay modes, which are reported in  The nonresonant contribution has been parametrized in an alternative way, P (m) = C(m)e a+bm . Each resonance parameter has been varied according to its uncertainty, and the meson radius used in the Blatt-Weisskopf damping factor [23], which is present in each relativistic Breit-Wigner term, has been varied between 0 and 2.5 GeV −1 . The amounts of backgrounds in the different final states have been varied according to their errors. The 3.9 GeV/c 2 structure in the DD mass spectrum has been alternatively described by a P-wave relativistic  Breit Wigner with free parameters. This effect dominates the systematic uncertainty on the ψ(4040) rate in the DD mass spectrum. The deviations from the central value are added in quadrature. Systematic effects also include the uncertainty on the total cross sections.
The corrected yields can also be used to compute the branching fraction ratios. The results are shown in Table VIII together with predictions of models: significant discrepancies are observed, expecially with the 3 P 0 model [5].
The D * D and D * D * mass spectra have been refit with an additional Y (4260) resonance, which is allowed to interfere with all the other terms.
The fit gives Y (4260) fractions of (2.2 ± 2.9 stat ± 2.5 syst )% and (4.0 ± 2.0 stat ± 4.2 syst )% corresponding to 18 ± 24 stat ± 21 syst and 9 ± 5 stat ± 10 syst events for Y (4260) → D * D and Y (4260) → D * D * , respectively. Systematic errors due to uncertainties on masses and widths of the ψ(4040), ψ(4160), ψ(4415), and Y (4260) resonances are evaluated by varying the masses and widths by their uncertainty in the fit. The amount of background in each final state is varied within its statistical error, and the meson radii in Breit-Wigner terms are varied between 0 and 2.5 GeV −1 . Deviations from the central value are added in quadrature.
These Y (4260) yields in the D * D and D * D * channels are used to compute the cross section times branching fraction, which can then be compared to our measure-ment from the J/ψπ + π − channel [2]. We obtain and at the 90% confidence level. Using the DD cross section measured in the earlier BABAR work [15], we obtain the sum of the e + e − → DD, e + e − → D * D, and e + e − → D * D * cross sections shown in Fig. 13: the arrow indicates the position of the Y (4260), which falls in a local minimum, in agreement with the cross section measured for hadron production in e + e − annihilation [9].

XI. CONCLUSIONS
We have studied the exclusive ISR production of the DD, D * D, and D * D * systems. The mass spectra show production of the J P C = 1 −− states ψ(3770), ψ(4040), ψ(4160), and ψ(4415). Fits to the mass spectra provide amplitudes and relative phases for the charmonium states, from which first measurements of branching fraction ratios are obtained. Finally, upper limits on Y (4260) → D * D and Y (4260) → D * D * decays are computed.
If the Y (4260) is a 1 −− charmonium state, it should decay predominantly to DD, D * D, and D * D * [5,6]. Within the present limited data sample size, no evidence is found for Y (4260) decays to DD, D * D, or D * D * . Other explanations for the Y (4260) have been proposed, such as a VIII: Ratios of branching fractions for the three ψ resonances. The first error is statistical, the second systematic. Theoretical expectations are from the 3 P0 model [5], C 3 model [6], and ρKρ model [14]. hybrid, baryonium, molecule or tetraquark state. In the case of a hybrid state, the decay rates to DD, D * D, and D * D * are expected to be small [10,24].

XII. ACKNOWLEDGEMENTS
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), the Commissariatà l'Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Education and Science of the Russian Federation, Ministerio de Educación y Ciencia (Spain), and the Science and Technology Facilities Council (United Kingdom). Individuals have received support from the Marie-Curie IEF program (European Union) and the A. P. Sloan Foundation.