Study of the Exclusive Initial-State-Radiation Production of the DDbar System

A search for charmonium and other new states is performed in a study of exclusive initial-state-radiation production of D Dbar events from electron-positron annihilations at a center-of-mass energy of 10.58 GeV. The data sample corresponds to an integrated luminosity of 384 fb-1 and was recorded by the BABAR experiment at the PEP-II storage ring. The D Dbar mass spectrum shows clear evidence of the psi(3770) plus other structures near 3.9, 4.1, and 4.4 GeV/c^2. No evidence for Y(4260) ->D Dbar is observed, leading to an upper limit of B(Y(4260) ->D Dbar)/B(Y(4260) ->J/psi pi+ pi-)<1.0 at 90 % confidence level.

PACS numbers: 13.66.Bc, 13.87.Fh,14.40.Gx The surprising discovery of new states decaying to J/ψπ + π − [1,2] has renewed interest in the field of charmonium spectroscopy, as the new states are not easy to accommodate in the quark model. In particular, the BABAR experiment has 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. Structure, possibly related to the Y (4260), has been observed in the reaction e + e − → γ ISR ψ(2S)π + π − [3]. A charmonium state at this mass would be expected to decay predominantly to DD, DD * or D * D * [4]. 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 [5], and that at the Y (4260) mass the cross section for e + e − → hadrons exhibits a local minimum [6]. Many theoretical interpretations for the Y (4260) have been proposed, including unconventional scenarios: quark-antiquark gluon hybrids [7], tetraquarks [8] and hadronic molecules [9]. For a discussion and a list of references see, for example, Ref. [10].
This work explores ISR production of the DD final state for evidence of charmonium states and unconventional structures. A study by the BELLE collaboration of the DD * , and D * D * final states can be found in Ref. [11].
This analysis is based on a 384 fb −1 data sample recorded at the Υ(4S) resonance and 40 MeV/c 2 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 [12]. Charged particles are detected and their momenta measured by 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. A ring-imaging Cherenkov detector (DIRC) combined with energy-loss measurements in the SVT and DCH are used to identify charged kaon and pion candidates. Photon energies are measured with a CsI(Tl) electromagnetic calorimeter (EMC).
DD candidates are reconstructed in seven combinations of D decay modes, listed in Table I [13]. In each channel we allow any number of photons in the event. Events are selected if the number of well-measured tracks is exactly equal to the total number of charged daughter particles for the D and the D final states. Neutral pion candidates are formed from pairs of photons each having an energy greater than 30 MeV. The K 0 S candidates are reconstructed in the π + π − decay mode. The tracks of each D candidate are geometrically constrained to come from a common vertex. Additionally, for the D 0 → K − π + π 0 channel, π 0 and D 0 mass constraints are included in the fit, and for the D − → K 0 S π − channel a K 0 S mass constraint is imposed. D candidates with a χ 2 fit probability greater than 0.1% are retained. Subsequently, each DD pair is refitted to a common vertex with the constraint that they originate from the e + e − interaction region; only candidates with a χ 2 fit probability greater than 0.1% are retained. Extra π 0 candidates may originate from random combinations of photons. Aside from π 0 's from D 0 decays, we require that there be no more than one other π 0 candidate in the event (except for channel 4, where we require that there are none).
For D channels without a π 0 candidate, the D momentum is determined from the summed 3-momenta of the decay particles and the energy is computed using the nominal D mass value [6,14]. For the D 0 → K − π + π 0 channel, the 4-momentum from the mass constrained fit is used. This procedure gives similar DD mass resolutions for all the channels.
The ISR photon, preferentially emitted at small angles with respect to the beam axis, escapes detection in approximately 90% of events. We therefore reconstruct the ISR photon as a missing particle. We define the squared recoil mass (M 2 rec ) to the DD system using the four-momenta of the beam particles (p e ± ) and the reconstructed D (p D ) and D (p D ): This quantity should peak near zero for ISR events and for exclusive production of e + e − → DD or e + e − → DD * . In the latter case, the DD mass distribution peaks at masses well above 6 GeV/c 2 . Therefore we select ISR events by requiring a DD invariant mass below 6 GeV/c 2 and | M 2 rec |< 1 GeV 2 /c 4 . Monte Carlo simulations of e + e − → γ ISR DD and candidates from the process e + e − → γ ISR J/ψ, J/ψ → K + K − π + π − in data are used to validate the requirement on the number of residual π 0 and the shape of the M 2 rec distribution.
To estimate the number of background events in the signal region, the two-dimensional space spanned by the invariant masses of the two D candidates in each event is divided into nine regions: a central signal region and eight sideband regions above and below the signal regions, as illustrated in Fig. 1 for DD candidates reconstructed for the case of the K − π + and K + π − modes. The mass range for the signal region is within ±2.5σ of the D mass, and the sideband regions are 2.5σ wide and are separated from the signal region by 3.5σ, where σ is the mass resolution determined from a fit of a single Gaussian to the D candidate mass spectrum.
The distribution of M 2 rec , summed over all DD channels, is shown in Fig. 2. The shaded histogram corresponds to the background in the signal region estimated from the DD mass sidebands. The small inset in Fig   shows the distribution of the DD center-of-mass polar angle θ for DD candidates with | M 2 rec |< 1 GeV 2 /c 4 . The sharp peak at cosθ = −1 is typical of ISR production and agrees with Monte Carlo simulations.
The purity of each reconstructed D channel is demonstrated in Fig. 3 where projections of the candidate D mass distribution for events with | M 2 rec |< 1 GeV 2 /c 4 are shown. Background is low in all channels.
The DD mass spectrum summed over all channels (860 events) is shown in Fig. 4 where the curves are the results from the fit described later. The shaded histogram represents the background determined using the DD sideband regions and corresponds to 17.5% and 7.1% of the signal candidates for D 0 D 0 and D + D − , respectively. We observe a clear ψ(3770) signal and other structures at the positions of ψ(4040) and ψ(4415). We also observe a significant structure in the 3.9 GeV/c 2 region, which may not be due to a resonance; the coupled-channel model of Ref. [15] in fact describes qualitatively the observed DD mass spectrum and the structure around 3.9 GeV/c 2 without any need for additional ψ states.
To understand the background, we compute the expected contribution from ISR production of the D * 0D0 system. Using Monte Carlo simulations and the cross section estimate from Ref. [11] we find ≈ 6 % as possible contamination. This is confirmed by the examination of Dγ and Dπ 0 mass distributions where we find little evidence for D * 0 signal. In contrast, strong evidence for D * 0 production is observed for M 2 rec > 1.5 GeV 2 /c 4 . We investigate the possibility of background contributions from DDX final states (where X = γ) by exploring events in the M 2 rec sideband region 1.5 < M 2 rec < 2.5 GeV 2 /c 4 . The DD mass spectrum for these events shows no structure. We conclude that the residual background rec |< 1 GeV 2 /c 4 and a DD invariant mass below 6 GeV/c 2 . (a) K − π + mass spectrum summed over channels 1, 2, and 3. (b) K − π + π 0 mass spectrum summed over channels 2 and 4. (c) K − π + π + π − mass spectrum summed over channels 3 and 4. (d) K − π + π + mass spectrum for channel 5.
to our signal is consistent with originating mostly from combinatorial non-DD events.
In order to measure efficiency and DD mass resolution, ISR events are simulated at eight different values of the DD invariant mass between 3.75 and 7.25 GeV/c 2 . These events are generated using the GEANT4 detector simulation package [16] and are processed through the same reconstruction and analysis chain as are real events. The mass-dependent efficiency for each channel is fitted using a second-order polynomial. The mass resolution is determined from the difference between generated and reconstructed DD mass. The DD mass resolution is similar for all channels and increases with DD mass from 1.5 to 5 MeV/c 2 . We observe good agreement between Monte Carlo and data M 2 rec distributions.
We define N i (m DD ) as the number of DD candidates for channel i. The channel branching fraction is B i , and ǫ i (m DD ) is the efficiency as parametrized by the fitted polynomial. We define as ǫ B i (m DD ) the product efficiency times branching fraction for each channel, and then compute ǫ B (m DD ) as The values of ǫ B i (m DD ) are proportional to the expected yield for each channel. Their values, integrated over the DD mass spectrum, are reported in Table I. The resulting yields, corrected for efficiency and branching fractions, are found to be consistent within the errors.
The DD cross section is computed using The differential luminosity is computed as [17] dL (5) where α is the fine-structure constant, x = 1 − m 2 DD /s, s is the square of the e + e − center-of-mass energy, m e is the electron mass, and L is the integrated luminosity of 384 fb −1 . The background-subtracted cross sections for D 0 D 0 and D + D − , averaged over 20 MeV/c 2 bins, are shown in Fig. 5. Systematic errors on the cross sections (10.9% for D 0 D 0 and 8.1% for D + D − ) include uncertainties in the particle identification efficiencies and tracking efficiency, possible inaccuracies in the simulation of extraneous π 0 candidates, and uncertainties in the background estimates (≈ 6 %) and on the luminosity function (≈ 1 %).
Integrating the cross sections in the ψ(3770) region (3.74-3.80 GeV/c 2 ), we compute the ratio of branching fractions, to be compared with the value of 1.28 ± 0.14 reported by the PDG [6]. We perform an unbinned maximum likelihood fit to the DD mass spectrum summed over all channels. The parameters of the ψ(4040), ψ(4160), and ψ(4415) are fixed to the values reported in Ref. [18] while the Y (4260) parameters are taken from our measurement from the J/ψπ + π − channel [2]. The parameters of the ψ(3770) are left free in the fit. In addition, we search for evidence of the Y (4260) in this spectrum. Resolutions effects have been ignored since the widths of the resonances are much larger than the experimental resolution.
We express the total DD production as where c i and φ i are free parameters, W i are spin-1 relativistic Breit-Wigner distributions, P represents the nonresonant contribution, B describes the non-DD background and f (0.829 ± 0.015) is the signal fraction. The efficiency ǫ B (m DD ) is almost linear and increases from ≈ 2 × 10 −3 to ≈ 4 × 10 −3 in the fitted mass region. It has been parametrized by a 2 nd order polynomial and it has been multiplied by P and W i . The data require that we include the 3.9 GeV/c 2 structure, as suggested in Ref. [15], which we parameterize empirically as the square root of a Gaussian times a phase factor ( √ Ge iφ2 ). The parameters of the Gaussian are left free, and the phase allows interference with the ψ states.
We find that, in order to have a satisfactory description of the data, interference must be allowed between the resonances and the non-resonant contribution P . The latter contribution is parametrized either as a linear (a + bm) or a threshold function (m − m th ) a e −bm−cm 2 , where m = m DD , m th is the threshold DD mass, and a, b and c are free parameters. This threshold function has also been used to describe the non-DD background B.
The two different parametrizations give similar results, which are considered in the evaluation of the systematic errors. These include also uncertainties in the D mass and on the overall DD mass scale. The size of the nonresonant production is determined by the fit.
The fit with a linear non-resonant contribution is shown in Fig. 4(a). Figure 4(b) shows an expanded view of the threshold region.
The systematic errors due to the masses and the widths of the ψ(4040), ψ(4160), ψ(4415) and Y (4260) resonances in the fit are evaluated by varying them by their statistical uncertainties. The signal fraction has been varied within its statistical error and the meson radius used in the Blatt-Weisskopf damping factor [19] present in the relativistic Breit-Wigner has been varied between 0 and 5 GeV −1 . The deviations from the central values are added in quadrature. The uncertainty on ǫ B (m DD ) is evaluated by using a weighted mean of branching fraction and efficiency uncertainties for the different channels. The fitted Y (4260) yield before efficiency correction is 0.2 ± 6.1 stat ± 2.8 syst events.
This Y (4260) yield in the DD channel is used to compute the cross section times branching fraction, which can then be compared to our measurement from the J/ψπ + π − channel [2]. We obtain or Γ(Y (4260) → e + e − ) · B(Y (4260) → DD) < 5.7 eV, (13) at 90% confidence level.
In conclusion, we have studied the exclusive ISR production of the DD system. The mass spectrum is dominated by J P C = 1 −− states; in particular the ψ(3770) is clearly seen. In order to fit the mass spectrum, signals from ψ(4040), ψ(4160), and ψ(4415) have been included. The fit requires the presence of a broad structure near 3900 MeV/c 2 . The presence of an enhancement in this region is predicted by a coupled channel model from Eichten et al. [15], although the possibility of the presence of a new ψ state cannot be excluded.
If the Y (4260) is a 1 −− charmonium state, it should decay predominantly to DD [4]; however no evidence is found for Y (4260) decays to DD. Other explanations have been proposed, such as a hybrid, baryonium, or tetraquark state.
We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE and NSF (USA), NSERC (Canada), CEA and CNRS-IN2P3 (France), BMBF and DFG (Germany), INFN (Italy), FOM (The Netherlands), NFR (Norway), MIST (Russia), MEC (Spain), and STFC (United Kingdom). Individuals have received support from the Marie Curie EIF (European Union) and the A. P. Sloan Foundation. * Deceased † Now at Tel Aviv University, Tel Aviv, 69978, Israel ‡ Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy § Also with Università della Basilicata, Potenza, Italy ¶ Also with Universitat de Barcelona, Facultat de Fisica, Departament ECM, E-08028 Barcelona, Spain