Evidence of a resonant structure in the $e^+e^-\to \pi^+D^0D^{*-}$ cross section between 4.05 and 4.60 GeV

The cross section of the process $e^+e^-\to \pi^+D^0D^{*-}$ for center-of-mass energies from 4.05 to 4.60~GeV is measured precisely using data samples collected with the BESIII detector operating at the BEPCII storage ring. Two enhancements are clearly visible in the cross section around 4.23 and 4.40~GeV. Using several models to describe the dressed cross section yields stable parameters for the first enhancement, which has a mass of $4228.6 \pm 4.1 \pm 6.3 \un{MeV}/c^2$ and a width of $77.0 \pm 6.8 \pm 6.3 \un{MeV}$, where the first uncertainties are statistical and the second ones are systematic. Our resonant mass is consistent with previous observations of the $Y(4220)$ state and the theoretical prediction of a $D\bar{D}_1(2420)$ molecule. This result is the first observation of $Y(4220)$ associated with an open-charm final state. Fits with three resonance functions with additional $Y(4260)$, $Y(4320)$, $Y(4360)$, $\psi(4415)$, or a new resonance, do not show significant contributions from either of these resonances. The second enhancement is not from a single known resonance. It could contain contributions from $\psi(4415)$ and other resonances, and a detailed amplitude analysis is required to better understand this enhancement.

The cross section of the process e + e − → π + D 0 D * − for center-of-mass energies from 4.05 to 4.60 GeV is measured precisely using data samples collected with the BESIII detector operating at the BEPCII storage ring. Two enhancements are clearly visible in the cross section around 4.23 and 4.40 GeV. Using several models to describe the dressed cross section yields stable parameters for the first enhancement, which has a mass of 4228.6±4.1±6.3 MeV/c 2 and a width of 77.0±6.8±6.3 MeV, where the first uncertainties are statistical and the second ones are systematic. Our resonant mass is consistent with previous observations of the Y (4220) state and the theoretical prediction of a DD1(2420) molecule. This result is the first observation of Y (4220) associated with an opencharm final state. Fits with three resonance functions with additional Y (4260), Y (4320), Y (4360), ψ(4415), or a new resonance, do not show significant contributions from either of these resonances. The second enhancement is not from a single known resonance. It could contain contributions from ψ(4415) and other resonances, and a detailed amplitude analysis is required to better understand this enhancement. As the first observed charmonium-like state with J PC = 1 −− , the Y (4260) has remained a mystery. Many experimental measurements and theoretical interpretations have been proposed for this state [1], such as hybrids [2], tetraquarks [3], and hadronic molecules [4].
Since it was observed only in hidden-charm processes, while its mass is close to open-charm thresholds, studies of the open-charm production cross section in e + e − annihilation will provide important information on its properties. The cross section for e + e − annihilation into D ( * )D( * ) pairs shows a dip at the resonance mass, 4.26 GeV/c 2 [5]. The Y (4260) mass is only about 29 MeV/c 2 below the nominal threshold for DD 1 (2420), which is the first open-charm relative S-wave channel coupling to J PC = 1 −− . The DD 1 (2420) molecule model is proposed as an interpretation of the Y (4260), but it predicts a significantly smaller mass of about 4.22 GeV/c 2 [6,7].
Recently, the precise measurement of the production cross section for e + e − → π + π − J/ψ from the BESIII experiment [8] indicates that the structure around 4260 MeV/c 2 actually consists of two resonances with masses 4222 MeV/c 2 and 4320 MeV/c 2 . The mass of the former resonance (referred to as Y (4220) hereafter) is consistent with the prediction of the DD 1 (2420) molecule model. Furthermore, a Y (4220) resonance has also been reported by the BESIII collaboration in the cross-section measurements of e + e − → ωχ c0 [9], e + e − → π + π − h c [10], and e + e − → π + π − ψ(3686) [11]. In addition, a new resonant structure with mass around 4.39 GeV/c 2 , the Y (4390), has been reported by BESIII in the reactions e + e − → π + π − h c [10] and e + e − → π + π − ψ(3686) [11]. The mass of the Y (4390) is about 45 MeV/c 2 and 70 MeV/c 2 higher than those of the Y (4360) [12] and the second component of the Y (4260) observed in e + e − → π + π − J/ψ by BESIII [8], respectively. The production of e + e − → πDD * is expected to be strongly enhanced above the nominal DD 1 (2420) threshold and could be a key for understanding existing puzzles with these Y states [7].
The cross section of e + e − → π + D 0 D * − was first measured by the Belle experiment using initial-state radiation (ISR) [13]. No evidence for charmonium(-like) states was found within their statistics. In this Letter, we report improved measurements of the production cross section of e + e − → π + D 0 D * − at center-of-mass energies ( √ s) from 4.05 to 4.60 GeV using data samples taken at 84 energy points [14] with the BESIII detector [15]. The data set contains 5 energy points ( √ s = 4.2263, 4.2580, 4.3583, 4.4156, and 4.5995 GeV) with integrated luminosities larger than 500 pb −1 ("H-XYZ data" hereafter), and 79 energy points with integrated luminosities smaller than 200 pb −1 . The D 0 meson is reconstructed in the D 0 → K − π + decay channel. The bachelor π + produced directly in the e + e − annihilation process is also reconstructed, while the D * − is not reconstructed directly, but is inferred from energy-momentum conservation. Charge-conjugate modes are implied, unless otherwise noted.
The BESIII detector is described in detail elsewhere [15]. A Monte Carlo (MC) simulation based on GEANT4 [16] includes the geometric description of the BESIII detector and its response. For each energy point, we generate MC samples of the signal process, e + e − → π + D 0 D * − , and the isospin partner process, e + e − → π + D − D * 0 , according to phase space (PHSP MC). The effect of ISR is simulated with KKMC [17] with a maximum energy for the ISR photon corresponding to the π + D 0 D * − mass threshold. Possible background contributions are estimated with KKMC-generated 'inclusive' MC samples with integrated luminosities comparable to the H-XYZ data, where the known decay modes are simulated with EVTGEN [18] using branching fractions taken from the PDG [12], and the remaining unknown decays are simulated with the LUNDCHARM model [19].
The charged tracks are reconstructed with standard selection requirements [20] and used to reconstruct D 0 meson candidates from K − π + track pairs. If there is more than one D 0 candidate in an event (∼0.3%), we choose the one whose invariant mass M (K − π + ) is closest to the world-average D 0 mass, m(D 0 ) [12]. The signal region is defined as To select the bachelor π + , at least one extra charged track, which is not used in the D 0 candidate and has charge opposite to that of the reconstructed K − , is chosen with the same selection criteria as described above.
The e + e − → D * D * background events are rejected by vetoing any D 0 π + candidates satisfying M (D 0 π + ) < 2.03 GeV/c 2 . After the above requirements, the presence of a D * − meson is inferred by the invariant mass recoiling against the D 0 π + system, RM (D 0 π + ). To improve the mass resolution, the corrected recoil mass, Fig. 1. If there is more than one bachelor π + in the event, the one whose RM cor (D 0 π + ) is closest to the world average D * − mass, m(D * − ), is selected. A study of the inclusive MC samples shows that only the isospin partner process e + e − → π + D − D * 0 (BKG1, hereafter) has an enhancement around the D * − mass region in the RM cor (D 0 π + ) distribution. The shape of this background process is different at each energy point and is taken from the MC simulation. The RM cor (D 0 π + ) distribution of the remaining background processes (BKG2, hereafter) does not peak and can be described by a firstorder polynomial function.
An unbinned maximum likelihood fit to the RM cor (D 0 π + ) distribution is performed to determine the signal yields. The signal shape is derived from the MC shape convolved with a Gaussian function. The background shape is parameterized as a sum of the shape from the PHSP MC sample for BKG1 and a first-order polynomial function for BKG2. We perform a simultaneous fit to the RM cor (D 0 π + ) distributions for all data samples to determine the yields of signal and background. The mean values of the Gaussian smearing function are constrained to be the same for all energy points. A center-of-mass energy-dependent width of the Gaussian function is obtained by fitting the widths of the five H-XYZ data samples with a first-order polynomial function, where these five widths are obtained by separate fits to the corresponding RM cor (D 0 π + ) distributions. The widths at √ s < 4.2263 GeV are fixed to that at √ s = 4.2263 GeV, since the fitted widths are close to zero. Figure 1 shows the fit result at √ s = 4.5995 GeV. The signal region is defined as where ∆M is the mean value of the Gaussian function obtained from the fit. A sideband region, used below, is defined as 1.91 < RM cor (D 0 π + ) < 1.95 GeV/c 2 .
The Born cross sections (σ Born ) and dressed cross sections (σ dress ) at the individual energy points are calcu- lated using where N obs is the signal yield, L is the integrated luminosity [21], 1 + δ is the ISR correction factor [22], 1 is the correction factor for vacuum polarization [23], B(D 0 → K − π + ) = (3.93 ± 0.04)% [12], and ǫ is the detection efficiency. Values of all above variables are given in the supplemental material [14]. Efficiencies at √ s = 4.2263, 4.2580, 4.3583, 4.4156 and 4.5995 GeV are calculated with MC simulated data samples [24] that are generated by the data-driven BODY3 generator based on EVTGEN [18], taking into account the influence of pos- [20,25] and highly-excited D states in the π + D 0 or π + D * − systems). Since the BODY3 generator requires a large selected sample obtained from events in the signal region after subtracting the background contribution (estimated with the events in the sideband region for BKG2 and MC simulation from BKG1), it is only used for the five energy points with high luminosity. Efficiencies at the other energy points are estimated with PHSP MC samples, with appropriate uncertainties included later. The obtained Born cross sections, which are consistent with and more precise than those of Belle [13], are summarized in the supplemental material [14].
The systematic uncertainties in the cross section measurements are listed in Table I. The uncertainty in luminosity is 1.0% at each energy point [26]. The uncertainty in B(D 0 → K − π + ) is 1.0% [12]. The uncertainty in the ISR correction factor is 3.0% [22]. The uncertainties associated with the detection efficiencies include tracking and PID efficiencies (1.0% per track), D 0 and D * − mass window requirements, and signal MC model. The uncertainties associated with the D 0 and D * − mass windows are estimated by repeating the analysis with an altered mass window requirement; the relative changes in the cross sections are taken as systematic uncertainties. The uncertainties associated with the BODY3 signal MC model consist of three parts: the choice of binning, and the BKG1 and BKG2 subtractions. The uncertainty associated with the choice of binning is estimated by repeating the simulation with an altered bin size. The uncertainty associated with the BKG1 subtraction is studied by replacing the PHSP MC sample with MC samples of processes including the intermediate states e + e − →D * 2 (2460) 0 D * 0 ,D * 2 (2460) 0 → π + D − and e + e − → D 1 (2460) + D − , D 1 (2460) + → π + D * 0 . For the BKG2 uncertainty, we replace the sideband events with the inclusive MC sample when subtracting the background. The maximum relative changes on the detection efficiency are taken as the corresponding uncertainties. The total signal MC model uncertainty is the sum in quadrature of these three contributions. To estimate the uncertainties of the signal MC model for the lowluminosity data, we estimate the detection efficiencies for the five energy points of large luminosity with the PHSP MC samples; the resultant largest difference with respect to the nominal efficiencies, 5.3%, is assigned as the corresponding uncertainty for the low-luminosity energy points. The uncertainties associated with signal shape, background shape, and fit range in the signal yield extraction are determined by changing the signal shape to the pure MC shape, by changing the background function from a linear polynomial function to a second-order one, and by changing the fit range, respectively. Due to limited statistics, fit results at the energy points with low luminosity suffer large statistical fluctuations in such refits; thus, the largest systematic uncertainties from the five large luminosity data samples are adopted. Assuming no significant correlations between sources, the total systematic uncertainty is obtained as the sum in quadrature.
The dressed cross section, which includes vacuum polarization effects, is shown in Fig. 2. Two enhancements around 4.23 and 4.40 GeV, denoted hereafter as R 1 and R 2 , are clearly visible. A maximum likelihood fit to the dressed cross section is performed to determine the parameters of these two enhancements. Since the measured cross sections have asymmetric uncertainties for the data with low statistics, the likelihood is described by an asymmetric Gaussian function as discussed in Ref. [10]. In the fit, the total amplitude is described by the coherent sum of a direct three-body phase-space term for e + e − → π + D 0 D * − and two relativistic Breit-Wigner (BW) functions, representing the resonant structures R 1 and R 2 : where the three-body phase-space factor P (m) [12] is modeled as a fixed fourth-order polynomial function.
The factor B j (m) = √ 12πΓ el j Γj m 2 −M 2 j +iMj Γj with j = 1 or 2 is the BW function for a vector state, where Γ el j = Γ e + e − B(π + D 0 D * − ) j is the product of the electronic partial width and the branching fraction to π + D 0 D * − . Parameters c, M j , Γ j , Γ el j , and φ j are the free parameters in the nominal fit. The beam energy spread (1.6 MeV) is considered by convolving with a Gaussian function whose width is 1.6 MeV. Only statistical uncertainties on the dressed cross sections are considered in the fit. There are four solutions with the same fit quality [27] and identical resonance parameters for R 1 and R 2 , but different c, Γ el j and φ j , as listed in Table II. We also fit the dressed cross section with the coherent sum of one resonance and a phase-space term; the change of the likelihood value, ∆(−2lnL), with respect to that of nominal fit including two resonances is 124.3. Taking into account the change in the number of degrees of freedom (4), the statistical significance of the two-resonance model over the one-resonance model is estimated to be 10.5σ.
Belle has observed the ψ(4415) in ψ(4415) → DD * 2 (2460) [28] which can also decay to the final state considered in this analysis. Considering the observations of other charmonium(-like) states, models fixing the mass and width of R 2 to those of Y (4260), Y (4320), Y (4360), or ψ(4415) are also investigated and ruled out with a confidence level equivalent to more than 5.0σ. Models including one additional known resonance, either Y (4260), Y (4320), Y (4360), or ψ(4415), in which the masses and widths of these resonances are fixed to the world average values [12], can improve the fit quality. However, the statistical significances of the additional resonances are only 0.4σ, 0.4σ, 1.4σ, and 1.0σ, respectively. The statistical significance of an additional unknown resonance is only 0.8σ, accounting for the two extra free parameters of mass and width. The high-mass enhancement has a more complicated underlying structure, the understanding of which requires a detailed amplitude analysis that is beyond the scope of this Letter.
All above models yield a stable set of parameters for R 1  but wildly varying parameters for R 2 , so we only estimate the systematic uncertainties of parameters of R 1 , which are mainly from the uncertainties of the absolute centerof-mass energy measurements, the cross section measurements, and the parameterization of the three-body phase-space factor. The uncertainty of the energy measurement (0.8 MeV) is propagated to the masses of the resonances. The uncertainty associated with the cross section measurements consists of two parts. The first is from the common uncertainties of the measured cross sections (tracking, PID, luminosity, and B(D 0 → K − π + )) for all energy points (4.5%); we shift up or down all measured cross sections by 4.5% simultaneously and repeat the fit on the measured cross sections. The differences, 0.1 MeV/c 2 for the mass and 0.1 MeV for the width, are taken as systematic uncertainties for the R 1 resonance. The second part includes all the other uncertainties of the measured cross sections. We add these uncertainties into the statistical ones in quadrature, and repeat the fit. The resulting differences in resonance parameters, 4.9 MeV/c 2 for the mass and 2.7 MeV for the width of R 1 , are taken as systematic uncertainties. The uncertainty associated with the three-body phase-space factor is determined by changing the parameterization function from a fourth-order polynomial function to a third-order one. The resulting differences, 3.8 MeV/c 2 for the mass and 5.7 MeV for the width, are taken as the corresponding uncertainties for R 1 . Assuming the individual systematic uncertainties are uncorrelated, the total systematic uncertainty is obtained by summing the individual values in quadrature, yielding 6.3 MeV/c 2 for the mass and 6.3 MeV for the width of R 1 .
In summary, the Born cross section for the process e + e − → π + D 0 D * − is precisely measured using the data samples collected at 84 energy points from 4.05 to 4.60 GeV with the BESIII detector. Two enhancements are observed in the dressed cross sections around 4.23 and 4.40 GeV. Using many models to describe the dressed cross section, we obtain a stable resonant structure around 4.23 GeV, the parameters of which are measured to be M (R 1 ) = (4228.6 ± 4.1 ± 6.3) MeV/c 2 and Γ(R 1 ) = (77.0 ± 6.8 ± 6.3) MeV, where the first uncertainties are statistical and the second ones systematic. The resonance parameters for the enhancement around 4.40 GeV are strongly dependent on the model assumptions, necessitating further studies.
The statistical significance of the two-resonance model over a one-resonance model is estimated to be 10.5σ. This is the first experimental evidence for open-charm production associated with the Y states. The mass of R 1 is consistent with the mass of the resonance observed in the hidden-charm processes by the BESIII experiment as well as the theoretical prediction of the DD 1 (2420) molecule interpretation [6]. The width of R 1 is consistent with that of e + e − → π + π − h c [10] and e + e − → π + π − ψ(3686) [11], but it is about 39 and 33 MeV/c 2 higher than that seen in e + e − → ωχ c0 [9] and e + e − → π + π − J/ψ [8], respectively. The minimum and maximum of the branching ratio, B(Y (4220)→πDD * ) B(Y (4220)→ππJ/ψ) ( B(Y (4220)→πDD * ) B(Y (4220)→ππhc) ), are calculated to be 1.3 ± 0.3 and 124.3 ± 36.1 (3.7 +2.5 −1.5 and 43.3 +29.0 −16.4 ) by assuming isospin symmetry, respectively. The measured Born cross section of e + e − → π + D 0 D * − at the Y (4220) peak is higher than the sum of the known hidden-charm channels. Since no other open-charm production associated with this Y state has yet been reported, the π + D 0 D * − final state may be the dominant decay mode of the Y (4220) state, as predicted by the DD 1 (2420) molecule interpretation [6]. No significant contributions from a third resonance are observed using three-resonance models with additional Y (4260), Y (4320), Y (4360), ψ(4415), or a new resonance, while Y (4320) and ψ(4415) are observed in e + e − → π + π − J/ψ [8] and ψ(4415) → DD * 2 (2460) [28], respectively. The amplitude studies of this final state and more studies on other open-charm production modes will shed additional light on the nature of these charmonium(-like) states.

Supplemental material
The integrated luminosities L, the numbers of signal events N obs , the ISR correction factors 1 + δ, the correction factors from vacuum polarization 1 |1−Π| 2 , the detection efficiencies ǫ, and the Born cross sections σ Born are summarized in Table III and Table IV. 15 energy points with integrated luminosities larger than 40 pb −1 are called "XY Z data" , and 69 energy points with integrated luminosities smaller than 20 pb −1 called "R-scan data" .