Observation of $e^{+}e^{-}\rightarrow \pi^{+}\pi^{-}\psi(3770)$ and $D_{1}(2420)^{0}\bar{D}^{0} + c.c.$

Several intermediate states of the reaction channels $e^{+}e^{-} \rightarrow \pi^{+} \pi^{-} D^{0} \bar{D}^{0}$ and $e^{+}e^{-} \rightarrow \pi^{+} \pi^{-} D^{+}D^{-}$ are studied using the data samples collected with the BESIII detector at center-of-mass energies above 4.08 GeV. For the first time in this final state, a $\psi(3770)$ signal is seen in the $D\bar{D}$ invariant mass spectrum, with a statistical significance of $5.2\sigma$ at $\sqrt{s} = 4.42$ GeV. There is also evidence for this resonance at $\sqrt{s}$ = 4.26 and 4.36 GeV with statistical significance of 3.2$\sigma$ and 3.3$\sigma$, respectively. In addition, the Born cross section of $e^{+}e^{-}\to \pi^{+}\pi^{-}\psi(3770)$ is measured. The proposed heavy-quark-spin-symmetry partner of the $X(3872)$, the state $X_{2}(4013)$, is also searched for in the $D\bar{D}$ invariant mass spectra. No obvious signal is found. The upper limit of the Born cross section of the process $e^{+}e^{-}\to \rho^{0}X_{2}(4013)$ combined with the branching fraction is measured. Also, the processes $e^{+}e^{-}\to D_{1}(2420)\bar{D} + c.c.$ are investigated. The neutral mode with $D_{1}(2420)^{0}\to D^{0}\pi^{+}\pi^{-}$ is reported with statistical significance of 7.4$\sigma$ at $\sqrt{s} = 4.42$ GeV for the first time, and evidence with statistical significance of 3.2$\sigma$ and 3.3$\sigma$ at $\sqrt{s}$ = 4.36 and 4.60 GeV is seen, respectively. No evident signal for the process $e^{+}e^{-}\to D_{1}(2420)^{0}\bar{D}^{0} + c.c., D_{1}(2420)^{0}\to D^{*+}\pi^{-}$ is reported. Evidence for $e^{+}e^{-}\to D_{1}(2420)^{+}D^{-} +~c.c., D_{1}(2420)^{+}\to D^{+}\pi^{+}\pi^{-}$ is reported with statistical significance of 3.1$\sigma$ and 3.0$\sigma$ at $\sqrt{s}$ = 4.36 and 4.42 GeV, respectively.

Several intermediate states of the reaction channels e + e − → π + π − D 0D0 and e + e − → π + π − D + D − are studied using the data samples collected with the BESIII detector at center-of-mass energies above 4.08 GeV. For the first time in this final state, a ψ(3770) signal is seen in the DD invariant mass spectrum, with a statistical significance of 5.2σ at √ s = 4.42 GeV. There is also evidence for this resonance at √ s = 4.26 and 4.36 GeV with statistical significance of 3.2σ and 3.3σ, respectively. In addition, the Born cross section of e + e − → π + π − ψ(3770) is measured. The proposed heavy-quark-spin-symmetry partner of the X(3872), the state X2(4013), is also searched for in the DD invariant mass spectra. No obvious signal is found. The upper limit of the Born cross section of the process e + e − → ρ 0 X2(4013) combined with the branching fraction is measured. Also, the processes e + e − → D1(2420)D + c.c. are investigated. The neutral mode with D1(2420) 0 → D 0 π + π − is reported with statistical significance of 7.4σ at √ s = 4.42 GeV for the first time, and evidence with statistical significance of 3.2σ and 3.3σ at √ s = 4.36 and 4.60 GeV is seen, respectively. No evident signal for the process e + e − → D1(2420) 0D0 + c.c., D1(2420) 0 → D * + π − is reported. Evidence for e + e − → D1(2420) + D − + c.c., D1(2420) + → D + π + π − is reported with statistical significance of 3.1σ and 3.0σ at √ s = 4. 36

I. INTRODUCTION
Heavy quarkonia have been studied for more than forty years for testing and developing quantum chromodynamics (QCD). On the one hand, some effective theories have been developed to describe quarkonium spectroscopy and transition dynamics [1][2][3]; on the other hand, many XY Z particles were discovered [4][5][6], and some of them are beyond the scope of potential models. The rich information gained from the XY Z particles may have opened a door through which quark confinement can be understood [7,8]. To understand these XY Z particles, it is of great importance to understand also the properties of the conventional quarkonia.
The ψ(3770) is generally assumed to be the 1 3 D 1 charmonium state with some admixture of the 2 3 S 1 state. One of the D-wave spin-triplet charmonium states, the ψ(1 3 D 2 ) or X(3823), has recently been observed in e + e − → π + π − ψ(1 3 D 2 ) at BESIII [17]. Therefore, the final states π + π − ψ(3770) and π + π − ψ(1 3 D 3 ) should be produced at BESIII as well, although so far there is no calculation on how large the production rates could be. The ψ(3770) decays dominantly to DD, which is also expected to be an important decay mode of the ψ(1 3 D 3 ). The predicted mass of the ψ(1 3 D 3 ) is at 3849 MeV/c 2 [18], however, there is no prediction for the width. Therefore, by studying the process e + e − → ππDD, one can also search for the ψ(1 3 D 3 ).
The X(3872) state was first observed by Belle [19], and confirmed subsequently by several other experiments [20][21][22]. Even though it clearly contains a cc pair, the X(3872) does not fit in the conventional charmonium spectrum. It could be interpreted as a DD * molecule with J P C = 1 ++ [23,24]. Throughout this paper, the charged conjugate mode is implied unless it is stated otherwise. Within this picture the existence of its heavy quark-spin-symmetry partner X 2 (4013) (J P C = 2 ++ ), an S-wave D * D * bound state, is predicted [25,26]. Its mass and width are predicted as about 4013 MeV/c 2 and ∼2-8 MeV, respectively. The X 2 (4013) is expected to decay dominantly to DD or DD * in D-wave. So it may also be produced in e + e − → π + π − DD. The possible discovery of the 2 ++ charmonium-like state will provide a strong support for the interpretation that the X(3872) is dominantly a DD * hadronic molecule [27].
Amongst various models to interpret the Y (4260) [9,11], the authors of Ref. [28] argue that the Y (4260) as a relative S-wave D 1 (2420)D system is able to accommodate nearly all the present observations of the Y (4260). Especially its absence in various open charm decay channels and the observation of the Z c (3900) in Y (4260) → π + π − J/ψ support this interpretation. In this model, the coupling strength of D 1 (2420)D to Y (4260) is a key piece of information. Because of D 1 (2420) decays to Dππ or D * π, this can also be studied via the ππDD final state.

II. THE EXPERIMENT AND DATA SETS
The BESIII detector is a magnetic spectrometer [29] located at the Beijing Electron Positron Collider (BEPCII) [30]. The cylindrical core of the BESIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
For this analysis, the data sets above 4.08 GeV recorded with the BESIII detector are used. The c.m. energy and the corresponding integrated luminosity of each data sample are listed in Table I. The c.m. energy is measured using di-muon events with a precision of 0.8 MeV [31]. The integrated luminosity is determined by analyzing large-angle Bhabha scattering events. The uncertainty of the integrated luminosity is 1.0% [32].
Simulated data samples produced with the GEANT4based [33] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the background contributions. The simulation includes the beam energy spread and allows for the production of initial state radiation (ISR) photons in the e + e − annihilation process. Both effects are modeled within the generator package KKMC [34].

III. EVENT SELECTION AND BACKGROUND ANALYSIS
In this analysis the DD (denoting D 0D0 and D + D − in the following) pairs are selected with both D mesons fully reconstructed in a number of hadronic decay channels (also called "double D tag" in the following). D 0 mesons are reconstructed in four decay modes (K − π + , K − π + π 0 , K − π + π + π − , and K − π + π + π − π 0 ) and D + mesons in five decay modes (K − π + π + , K − π + π + π 0 , K 0 S π + , K 0 S π + π 0 , and K 0 S π + π − π + ). TheD 0 and D − mesons are reconstructed in the charge conjugate final states of the D 0 and D + mesons, respectively. One π + π − pair is selected in addition to the tracks from DD decays.
Charged tracks are reconstructed from MDC hits within a polar-angle (θ) acceptance range of | cos θ| < 0.93 and required to pass within 10 cm of the interaction point in the beam direction and within 1 cm in the plane perpendicular to the beam. The TOF and dE/dx information is combined for each charged track to calculate the particle identification (PID) probability P i (i = π, K) of each particle-type hypothesis. P K > P π is required for a kaon candidate and P π > P K is required for a pion candidate. Tracks used in reconstructing K 0 S are exempted from these requirements. Electromagnetic showers are reconstructed by clustering the energy deposits of the EMC crystals. Efficiency and energy resolution are improved by including energy deposits in nearby TOF counters. A photon candidate is defined as a shower with an energy deposit of at least 25 MeV in the "barrel" region (| cos θ| < 0.8), or at least 50 MeV in the "endcap" region (0.86 < | cos θ| < 0.92). Showers in the transition region between the barrel and the end-cap are not well measured and are rejected. An additional requirement on the EMC hit timing (0 ≤ T ≤ 700 ns relative to the event start time) suppresses electronic noise and energy deposits unrelated to the event. To eliminate showers from bremsstrahlung photons which originated from charged particles, the angle between the shower and nearest charged track is required to be greater than 20 degrees. π 0 candidates are reconstructed from pairs of photon candidates with an invariant mass in the range 0.115 < M γγ < 0.150 GeV/c 2 . A one-constraint (1C) kinematic fit with the mass of the π 0 constrained to the world average value [37] is performed to improve the energy resolution.
K 0 S candidates are reconstructed from two oppositely charged tracks which satisfy | cos θ| < 0.93 for the polar angle and the distance to the average beam position in beam direction within 20 cm. For each pair of tracks, assuming they are π + and π − , a vertex fit is performed and the resulting track parameters are used to obtain the ππ invariant mass which must be in the range 0.487 < M ππ < 0.511 GeV/c 2 . The χ 2 from the vertex fit is required to be smaller than 100.
The selected K ± , π ± , K 0 S , and π 0 candidates are used to reconstruct D meson candidates which are composed to D 0D0 and D + D − meson pairs. If more than one DD pair per event is found with both D mesons decaying in the same way, the pair with the average massM = [M (D) + M (D)]/2 closest to the nominal mass of the D meson [37] is chosen. In each event, one negative and one positive charged π are required in addition. To reduce the background contribution and improve the mass resolution, a four-constraint (4C) kinematic fit is performed. The total four-momentum of all selected charged tracks and good photons from π 0 are constrained to that of the initial e + e − system. If the final state contains a π 0 or K 0 S meson, its mass is constrained in the kinematic fit as well. If there are multiple candidates in an event, the one with the smallest χ 2 of the kinematic fit is chosen. To find the optimal χ 2 criteria, the figure of merit FOM = ns √ ns+n b is maximized. Here n s is the number of signal events from signal MC simulation and n b is the number of backgrounds events from inclusive MC samples. The χ 2 is required to be less than 56 for the π + π − D 0D0 final state with selection efficiency of 90.1% and background rejection rate of 45.5%, and less than 40 for the π + π − D + D − final state with selection efficiency of 90.3% and background rejection rate of 29.5%.
The inclusive MC sample is used to investigate possible background contributions. There is neither a peaking background contribution found near 3.773 GeV/c 2 and 4.013 GeV/c 2 in the DD invariant mass distribution nor near 2.42 GeV/c 2 in the Dππ invariant mass distribution. From the study of the MC samples with highly excited charmed mesons, we find that only the process e + e − → D * 2 (2460) 0D0 , D * 2 (2460) 0 → D * + π − produces a peak near 2.46 GeV/c 2 in the signal region of the D 1 (2420) 0 in the invariant mass distribution of D * + π − . Therefore, There will be some non-DD backgrounds remaining in the signal region. According to the study of the inclusive MC, in the DD and Dππ invariant mass distribution, non-DD backgrounds and sidebands events are consistent with each other. Therefore, the events from the sidebands are used to describe non-DD backgrounds in this analysis.

IV. SIGNAL YIELD DETERMINATION
After imposing all the requirements mentioned above, the DD invariant mass distributions are shown in Fig To determine the signal yield of e + e − → π + π − ψ(3770), an unbinned maximum likelihood fit is performed to the M (DD) spectra as shown in Fig. 3. The signal contribution is described by the MC simulated shape which is modeled using non-parametric kernel-estimation [40]. The background component includes the channel D 1 (2420)D, the four-body decay π + π − DD (both described with the shape taken from the MC simulation which are also modeled with non-parametric kernel-estimation) and the non-DD background distribution (described by the DD sideband events). In the fit, the signal yields are free parameters with lower limits set to 0. The yields of the D 1 (2420)D and the π + π − DD background contributions are free parameters. The number of non-DD background events is fixed to the number of events from the sidebands. The signal yields at √ s = 4.26, 4.36, and 4.42 GeV returned by the fit are 30.7 ± 9.9, 68.7 ± 21.8, and 99.2 ± 21.0 events, respectively. The statistical significance of the signal yield is determined to be 3.2σ, 3.3σ, and 5.2σ, respectively, by comparing the log-likelihood values with and without the 1 and the green long-dashed lines the distributions from π + π − ψ(3770). The pink dotted-dotted-dashed lines show the distribution from π + π − X2(4013). The distributions from D 0 D 0 1 are normalized to the maximum bin content of data, and the distributions from π + π − ψ(3770) and π + π − X2(4013) are normalized arbitrarily. signal hypothesis and taking the change in the number of degrees of freedom into account. With the same method, the data samples taken at other c.m. energies are also studied as shown in Fig. 8 of Appendix A. The signal yields are listed in Table I. In this analysis, if the statistical significance of the signal is less than 1σ, we will scan the likelihood distribution as a function of the signal yield greater than 0 and use the difference between the most probable values and the thresholds of 68.3% total integral area as the errors.
We also search for structures in the π ± ψ(3770) invariant mass distribution at the energy points where the ππψ(3770) signal is most prominent. The π ± ψ(3770) distribution after requiring M (DD) ∈ [3.73, 3.80] GeV/c 2 around the ψ(3770) mass is shown in Fig. 4. There are hints for peaks at 4.04 and 4.13 GeV/c 2 in √ s = 4.42 GeV data, but the statistical significance is low.
For the search for the X 2 (4013) resonance, the region of large DD invariant masses is investigated. Figure 5 shows the distributions after imposing all the requirements above. There is no obvious signal visible around 4.013 GeV/c 2 . We try to fit these distributions with the signal shape of the process ρ 0 X 2 (4013) taken from the MC simulation and a third order polynomial as background distribution as shown in Fig. 5 Table II. C. e + e − → D1(2420)D After imposing all the requirements above, the Dπ + π − invariant mass distributions are shown in Fig. 6. A peak at  To determine the signal yield for the reaction e + e − → D 1 (2420) 0D0 , the D 0 π + π − invariant mass distribution is fitted with the signal shape taken from the MC simulation convolved with a Gaussian function to take into account the shift of the reconstructed mass to the generated one and the difference in the mass resolution between data and MC simulation as shown in Fig. 6. As background components, the channels π + π − ψ(3770) and π + π − D 0D0 are included in the fit as well as a non-DD component, which is fixed to the shape and number of events expected from the sideband regions. The numbers of π + π − ψ(3770) events are fixed to the values calculated using the cross sections we measured (see Table II). The yields of the signal events and of the π + π − D 0D0 events are allowed to vary in the fit. A similar fit is performed to the D * + π − invariant mass distribution as shown in Fig. 6. The signal shape is taken from the MC simulation in the same way as described above. As background components, the channels D * D π and D * 2 (2460) 0D are included in the fit as well as the non-DD component taken from the sideband regions. The signal yields at √ s = 4.36, 4.42, and 4.60 GeV are 17.8 ± 9.3, 22.3 ± 13.2, and 12.6 ± 7.3 events with the statistical significance of 1.6σ, 2.4σ, and 1.5σ, respectively.
To determine the signal yield of e + e − → D 1 (2420) + D − , the D + π + π − invariant mass distribution is fitted with a procedure similar to the neutral mode as shown in Fig. 6 The data samples taken at other c.m. energies are also studied with the same method. The fits are shown in Figs. 9-11 in Appendix B. Signal yields are listed in Tables III-V.

V. CROSS SECTION RESULTS
The Born cross section of e + e − → π + π − ψ(3770) is calculated with and B C are the branching fractions for ψ(3770) → D 0D0 and ψ(3770) → D + D − and B i (B j ) is the branching fraction for D 0 → i (D 0 → j) taken from PDG [37]. The same applies to ǫ k,l and B k (B l ) for charged mode. f v = 1 |1−Π| 2 is the vacuum polarization factor [41] and f r = (1 + δ r ) is the radiative correction factor which is defined as F (x, s) is the radiator function, which is calculated from QED [42] with an accuracy of 0.1%, γ is the ISR photon energy and m is the invariant mass of the final state after radiating the photon. σ dressed (s) is the energy dependent dressed cross section of e + e − → π + π − ψ(3770). Here the observed signal events are assumed to originate from the Y (4260) resonance to obtain the efficiency and the ISR correction factor. Then the mea-sured line shape is used as input to calculate the efficiency and ISR correction factor again. This procedure is repeated until the difference between two subsequent iterations is comparable with the statistical uncertainties. The Born cross sections are listed in Table I and shown in Fig. 7. At the energy points where no significant ψ(3770) signal yields are observed (significance < 5σ) the upper limits on the cross sections at the 90% confidence level (C.L.) are calculated using the Bayesian method [37] with a flat prior. By fitting the DD invariant mass distribution with fixed values for the signal yield, we obtain a scan of the likelihood distribution as a function of the cross section. To take all systematic uncertainties into consideration we convolve the likelihood distribution with a Gaussian function with a width corresponding to the total systematic uncertainty. The upper limit on σ at the 90% C.L. is obtained TABLE I: Results for the process e + e − → π + π − ψ(3770). Shown in the table are the number of observed events N obs , the integrated luminosity Lint, the radiation correction factor 1 + δ r , the vacuum polarization correction factor 1 |1−Π| 2 , the summation over the products of branching fraction and efficiency i,j ǫi,jBiBj (left) for the D 0D0 and (right) for the D + D − mode, the Born cross section σ B and the statistical significance S. The upper limits are at 90% C.L.     TABLE III: Results for the process e + e − → D1(2420) 0D0 with D1(2420) 0 → D 0 π + π − + c.c. For the symbols see Table I, the penultimate column is the Born cross section σ B times the branching fraction of D1(2420) 0 → D 0 π + π − . The upper limits are at 90% C.L.
All upper limits on the cross sections are also listed in Table I.  Results for the reaction channel e + e − → D1(2420) 0D0 with D1(2420) 0 → D * + π − + c.c. (For symbols see Table III).  Results for the reaction channel e + e − → D1(2420) + D − with D1(2420) + → D + π + π − + c.c. (For symbols see Table III). 90% C.L. are estimated using the same method as described above. All results and upper limits are listed in Table II. For the reaction channel e + e − → D 1 (2420)D with D 1 (2420) → X(Dπ + π − or D * π), the product of the Born cross section times the branching fraction of D 1 (2420) → X is calculated using where ǫ i,j is the selection efficiency for each process e + e − → D 1 (2420)D (D 1 (2420) → Dππ/D * π, D → i,D → j). The low momentum of the π meson from D * → Dπ decay reduces the efficiency for the decay channel D 1 (2420) → D * π in comparison to D 1 (2420) → Dππ. The other variables are the same as defined in Eq. (1). For the D 1 (2420) → D * π channel, the branching fraction B D * →Dπ is taken into account. The cross sections as a function of c.m. energy are shown in Fig. 7. At the energy points where no significant D 1 signals are observed (significance < 5σ), the upper limits on the cross sections at the 90% C.L. are estimated using the same method as described above. All numbers are shown in the Tables III, IV, and V for the neutral and the charged modes, respectively.

VI. SYSTEMATIC UNCERTAINTY ESTIMATION
The systematic uncertainties in the cross section measurements mainly stem from the integrated luminosity, the tracking and photon detection efficiency, various intermediate branching fractions, the ISR correction factor, the signal and background shapes, the fit range, the signal region of double tag of the D mesons, and the cross section of the π + π − ψ(3770) final state. The estimation of the systematic uncertainties is described in the following and the results at  (a) The uncertainty from the integrated luminosity measurement using Bhabha (e + e − → e + e − ) scattering events is esti-   mated to be 1.0% [32].
(c) ISR photons are simulated by using the KKMC package. The shape of the c.m. energy dependent cross section affects the radiative correction factor and the reconstruction efficiency. For the reactions e + e − → π + π − ψ(3770) and e + e − → D 1 (2420)D, the difference between the last two iterations is taken as the systematic uncertainty. Since we have no knowledge on the production cross section for the reaction e + e − → ρ 0 X 2 (4013), we assume that the cross section of ρ 0 X 2 (4013) follows the Y (4260) or the ψ(4415) line shape. The difference between these two assumptions is taken as the systematic uncertainty.
(d) For the determination of the systematic uncertainty caused by the signal shape, additional MC samples are produced by varying the width of the signal resonance by one standard deviation of its world average value [37]. The largest difference of the cross section compared with the nominal value is taken as the systematic uncertainty of the signal shape.
(e) The systematic uncertainty caused by the background shape, which is taken from MC simulation of the final states π + π − ψ(3770), D 1 (2420)D and D 1 (2460)D, is estimated by generating alternative MC samples where the width of the ψ(3770), D 1 (2420) and D 1 (2460) resonances is changed by one standard deviation of the world average value [37]. The largest difference of the cross section compared with the nominal fit value is taken as the systematic uncertainty of the background shape. The systematic uncertainty originating from the sideband selection is estimated by changing the sideband windows by 10 MeV/c 2 . The largest difference of the cross section compared with the nominal mass window is taken as systematic uncertainty. For e + e − → ρ 0 X 2 (4013), the background shape is changed from a third order polynomial to a fourth order polynomial, and the difference is taken as the systematic uncertainty of background shape.
(f) The systematic uncertainty caused by the choice of the fit range is estimated by varying the limits of the fit range by 20 MeV/c 2 . The largest difference of the cross section from the nominal value is taken as systematic uncertainty.
(g) In order to estimate the systematic uncertainty due to the selection of the signal window for the double D-tag method, the whole analysis is repeated by changing the signal region from |∆M | < 35 MeV/c 2 , −6 < ∆M < 10 MeV/c 2 to |∆M | < 39 MeV/c 2 , −8 < ∆M < 12 MeV/c 2 for the D 0D0 mode and from |∆M | < 25 MeV/c 2 , −5 < ∆M < 10 MeV/c 2 to |∆M | < 29 MeV/c 2 , −7 < ∆M < 12 MeV/c 2 for the D + D − mode. The difference of the cross section from the nominal value is taken as systematic uncertainty.
(h) The systematic uncertainty caused by the fixed number of π + π − ψ(3770) events in the fit of M (Dπ + π − ) is estimated by varying the fixed number by one standard deviation. The largest deviation from the nominal cross section is taken as systematic uncertainty.
Tables VI to X summarize all the systematic uncertainties. The overall systematic uncertainty for each process and c.m. energy is obtained by summing up all sources of systematic uncertainties in quadrature, assuming they are uncorrelated.
Whether the events are from the Y (4390) or the ψ(4415) resonance or from any other resonance cannot be distinguished based on the current statistics. For the data points with enough statistics for the e + e − → π + π − ψ(3770) final state, no significant structure (i.e. a possible Z c state) is observed in the π ± ψ(3770) system.