Measurement of $e^{+}e^{-} \to \pi^{0}\pi^{0}\psi(3686)$ at $\sqrt{s}$ from 4.009 to 4.600 GeV and observation of a neutral charmoniumlike structure

Using $e^{+}e^{-}$ collision data collected with the BESIII detector at the BEPCII collider corresponding to an integrated luminosity of 5.2 $fb^{-1}$ at center-of-mass energies ($\sqrt{s}$) from 4.009 to 4.600 GeV, the process $e^{+}e^{-} \to \pi^{0} \pi^{0}\psi(3686)$ is studied for the first time. The corresponding Born cross sections are measured and found to be half of those of the reaction $e^{+}e^{-} \to \pi^{+}\pi^{-}\psi(3686)$. This is consistent with the expectation from isospin symmetry. Furthermore, the Dalitz plots for $\pi^{0}\pi^{0}\psi(3686)$ are accordant with those of $\pi^{+}\pi^{-}\psi(3686)$ at all energy points, and a neutral analogue to the structure in $\pi^\pm \psi(3686)$ around 4040 MeV/$c^{2}$ first observed at $\sqrt{s}$=4.416 GeV is observed in the isospin neutral mode at the same energy.

process e + e − → π 0 π 0 ψ(3686) is studied for the first time. The corresponding Born cross sections are measured and found to be half of those of the reaction e + e − → π + π − ψ(3686). This is consistent with the expectation from isospin symmetry. Furthermore, the Dalitz plots for π 0 π 0 ψ(3686) are accordant with those of π + π − ψ(3686) at all energy points, and a neutral analogue to the structure in π ± ψ(3686) around 4040 MeV/c 2 first observed at √ s = 4.416 GeV is observed in the isospin neutral mode at the same energy. The vector charmoniumlike state Y (4360) was observed and subsequently confirmed in e + e − → (γ ISR )π + π − ψ(3686) by BaBar, Belle, and BESIII [1][2][3], where γ ISR refers to an Initial State Radiation (ISR) photon. However, the nature of the Y (4360) remains mysterious [4], as is the case for other states of the Y family, e.g. the Y (4260) observed in e + e − → (γ ISR )π + π − J/ψ [5][6][7][8]. Many theoretical interpretations have been proposed to explain the underlying structure of the Y family of states [9][10][11]. It is therefore compelling to study the Y (4360) in its π 0 π 0 transition to ψ(3686) and to examine isospin symmetry.
In recent years, a new pattern of charmoniumlike states, the Z ± c 's, was observed in the systems of a charged pion and a low-mass charmonium state [3,7,[12][13][14], as well as in charmed mesons pairs [15][16][17]. The observation of Z ± c particles and of similar states in the bottomonium system [18] indicates the discovery of a new class of hadrons [19]. More recently, neutral charmoniumlike states, which are referred to as Z 0 c 's, have been reported in analogous systems [20][21][22][23]. These are regarded as the neutral isospin partners of the Z ± c 's. A charmoniumlike structure observed in e + e − → π + π − ψ(3686) by BESIII [3] was also reported in Belle's latest updated result [2]. By analogy, it is interesting to search for its neutral isospin partner in e + e − → π 0 π 0 ψ(3686).
In this Letter, we present a study of the process e + e − → π 0 π 0 ψ(3686) at center-of-mass (c.m.) energies ( √ s) from 4.009 to 4.600 GeV. The corresponding Born cross sections are measured for the first time. A new neutral structure is observed in the π 0 ψ(3686) invariant mass spectra around 4040 MeV/c 2 . The data samples used in this analysis were collected with the BESIII detector at 16 different c.m. energies with a total integrated luminosity of 5.2 fb −1 [24]. The c.m. energies have been measured with di-muon events for each energy point [25].
The signal candidates are required to have four charged tracks with zero net charge and at least four photon candidates. The selection criteria for good charged tracks and photons, the separation between pions, electrons and muons, as well as the hit number required in the muon system for the µ + µ − pair are the same as those in Ref. [3].
A four-constraint (4C) kinematic fit imposing energymomentum conservation under the hypothesis e + e − → γγγγπ + π − ℓ + ℓ − is performed, and χ 2 4C < 120 is required. For the events with more than four photons, the combination with the least χ 2 4C is retained. The pairing of photons into the two π 0 is chosen by minimizing (M (γ 1 γ 2 ) − M (π 0 )) 2 + (M (γ 3 γ 4 ) − M (π 0 )) 2 . The J/ψ and π 0 candidates are selected by requiring 3.05 < M (ℓ + ℓ − ) < 3.15 GeV/c 2 and |M (γ i γ j ) − M (π 0 )| < 20 MeV/c 2 , where M (π 0 ) is the π 0 mass according to the PDG [34]. A seven-constraint (7C) kinematic fit with additional constraints on the two π 0 and J/ψ masses [34] is imposed to suppress the non-π 0 π 0 π + π − J/ψ backgrounds and improve the mass resolution. Figure 1 shows the scatter plots of the π + π − recoil mass M recoil (π + π − ) versus M (π + π − J/ψ) and the M (π + π − J/ψ) spectra for the data samples at √ s = 4.226, 4.258, 4.358, 4.416, and 4.600 GeV, which have relatively large statistics. Vertical and horizontal bands at the ψ(3686) mass position are observed clearly in the scatter plots, corresponding to the processes e + e − → π 0 π 0 ψ(3686) and e + e − → π + π − ψ(3686), respectively. The narrow peaks in the M (π + π − J/ψ) spectra indicate the signal process e + e − → π 0 π 0 ψ(3686), while the relative broad bumps with position depending upon c.m. energy are from e + e − → π + π − ψ(3686). The inclusive and exclusive MC samples, as well as the data in the J/ψ sideband region (selected by applying a 6C kinematic fit without the J/ψ mass constraint instead of the 7C kinematic fit), are used to investigate the backgrounds. The dominant background is e + e − → π + π − ψ(3686), which has the same final states as the signal. An unbinned maximum likelihood fit is performed to the M (π + π − J/ψ) spectra to determine the signal yields. In the fit, the probability density function (PDF) of e + e − → π 0 π 0 ψ(3686) and e + e − → π + π − ψ(3686) are described with the MC simulated shapes convolved with a Gaussian function, where the parameters of the Gaussian function are determined in the fit, in order to account for the resolution difference and potential mass shift between the data and MC simulation. The other backgrounds are described with a linear function. The fits curves are shown in Fig. 1.
The Born cross section is calculated from where L int is the integrated luminosity, N obs is the signal yield from the fit, (1 + δ r ) is the ISR correction factor which is obtained by using a QED calculation [32] and incorporating the input lineshape of the cross section which is taken to be the same as that of e + e − → π + π − ψ(3686) from the Belle experiment [2], (1 + δ v ) is the vacuum polarization factor taken from a QED calculation with an accuracy of 0.5% [33], B is the product of the branching fractions in the decay chain, 3.89%, taken from the PDG [34], and ǫ is the detection efficiency. The numbers used in the Born cross section calculation and the cross sections are summarized in the Supplemental Material [35]. The comparison of the Born cross section of e + e − → π 0 π 0 ψ(3686) to that of e + e − → π + π − ψ(3686) for the data samples with large luminosities is shown in Fig. 2. The Born cross sections of e + e − → π + π − ψ(3686) are also calculated with the corresponding event yields, and are consistent with the results in Ref. [3]. An alternative fit to the M recoil (π + π − ) spectra, which have a narrow peak for e + e − → π + π − ψ(3686) and a broad bump depending on c.m. energy for e + e − → π 0 π 0 ψ(3686), is performed. In the fit, the PDF is described by the similar strategy with the M (π + π − J/ψ) spectra. The resulting Born cross sections are consistent with the nominal values for the two processes.
Upper limits at the 90% confidence level (C.L.) on the Born cross sections are determined by using a frequentist method with profile likelihood treatment of systematic uncertainties [36]. The number of signal events (N obs ) is counted in the region 3.671 < M (π + π − J/ψ) < 3.701 GeV/c 2 , while the number of background events (N bkg ) is evaluated in the region 3.630 < M (π + π − J/ψ) < 3.660 GeV/c 2 or 3.712 < M (π + π − J/ψ) < 3.742 GeV/c 2 . In the calculation, the observed events are assumed to have a Poisson distribution and the event selection efficiencies to follow a Gaussian distribution. The upper limits are shown in the supplemental material [35]. The cross section ratios, R π + π − ψ(3686) = σ(e + e − →π 0 π 0 ψ(3686)) σ(e + e − →π + π − ψ(3686)) , are calculated for data samples with large luminosities and are listed in the Supplemental Material [35], where σ(e + e − → π + π − ψ(3686)) are taken from Ref. [3]. A set of common systematic uncertainties among the two processes, including those on luminosity, tracking efficiencies, and the requirements on the lepton tracks, cancel in the calculation. The weighted average of the ratios at √ s = 4.226, 4.258, 4.358, 4.416 GeV is 0.48 ± 0.04 ± 0.02. Within uncertainties, the resulting R π + π − ψ(3686) is consistent with the value 0.5 expected from isospin symmetry.
The following sources of systematic uncertainty are considered in the cross section measurements. The uncertainty on the efficiency for charged tracks (photons) is 1% per track (photon) [37,38]. The uncertainty on the hit number requirement in the muon counter is 4.2%, obtained by studying a sample of e + e − → π + π − J/ψ events. The uncertainty related with the kinematic fit is estimated with the same method as in Ref. [39], and is in the range 0.3% to 1.2% depending on c.m. energy. The uncertainties of the π 0 and J/ψ invariant mass requirements are evaluated by tuning the corresponding MC distributions according to data, and are in the range 0.2% to 0.5% and 0.1% to 0.3%, respectively, depending on c.m. energy. The uncertainties related to the fit procedure are investigated by varying the fit range, replacing the linear function for the background by a second-order polynomial function for background, and varying the width of the Gaussian function for the signal, and are in the range 1.6% to 7.3% depending on c.m. energy. For the data samples with large luminosity, the detection efficiencies are estimated with the MC samples re-weighted according to the Dalitz plots distributions of M 2 (π 0 ψ(3686)) versus M 2 (π 0 π 0 ) found in data. The corresponding uncertainty is estimated by varying the weighting factors according to the statistical uncertainty in each bin. For the data samples with low luminosity, the detection efficiencies are estimated with the Jpipi model MC samples. The corresponding systematic uncertainties are estimated with the data samples with large luminosity by comparing the efficiencies derived from the Jpipi model MC sample with the nominal model. The uncertainty associated with the ISR correction factor is studied by replacing the input cross section line shape with the latest results from BaBar [1] in the KKMC generator, and is in the range 0.3% to 2.4% depending on c.m. energy. The uncertainty of the vacuum polarization factor is 0.5% from a QCD approach [33].
The uncertainty of the integrated luminosity is 1%, as determined with large-angle Bhabha events [24]. The uncertainties of the branching fractions of the intermediate states are taken from the PDG [34]. A summary of all considered systematic uncertainties is shown in the supplemental material [35]. Assuming all sources of systematic errors are independent, the total uncertainties are the quadratic sums of the individual values, ranging from 7.8% to 10.8%, depending on the c.m. energy.
Possible intermediate states in e + e − → π 0 π 0 ψ(3686) are investigated using the data samples at √ s = 4.226, 4.258, 4.358, and 4.416 GeV. The ψ(3686) signal is extracted by selecting the events in the mass range 3.676 < M (π + π − ℓ + ℓ − ) < 3.696 GeV/c 2 . The Dalitz plots M 2 (π 0 π 0 ) versus M 2 (π 0 ψ(3686)) as well as the corresponding one-dimensional distributions are shown in Fig. 3. Good agreement of these distributions with those observed in the charged mode in Ref. [3] is found, which confirms the variations of the kinematic behavior at different energy points and demonstrates isospin conservation. A structure with a mass around 4040 MeV/c 2 in the M (π 0 ψ(3686)) spectrum is observed in the data sample at √ s = 4.416 GeV, while two bumps around 3900 and 4040 MeV/c 2 are evident in the data sample at √ s = 4.258 GeV. It is worth noting that for the data sample at √ s = 4.226 GeV, this structure is not observed in the M (π 0 ψ(3686)) distribution. The behavior observed is similar with that in the charged mode [3]. The dominant background is e + e − → π + π − ψ(3686) as shown in Fig. 1. The other backgrounds are found to be negligible from the study of sideband region.
An unbinned maximum likelihood fit is performed to the Dalitz plot of M 2 (π 0 1 ψ(3686)) versus M 2 (π 0 2 ψ(3686)) (denoted as x and y in Eq. (2), respectively) to determine the properties of the observed structure at √ s = 4.416 GeV. In the fit, the observed structure is assumed to be a neutral charmoniumlike state with spin-parity 1 + , modeled with an S-wave Breit-Wigner function in two dimensions, taking into account the mass resolution and detection efficiency, where, p 1/2 (q 1/2 ) is the momentum of the charmoniumlike state (ψ(3686)) in the rest frame of its mother particle, and M R and Γ are the mass and width of the charmoniumlike state, respectively. The PDF of the process e + e − → π 0 π 0 ψ(3686) without an intermediate state is taken from the Jpipi model MC simulation. The background is found to be negligible, and is not included in the fit. Since the two π 0 mesons in the final state are experimentally indistinguishable, the fit is performed with two entries per event, and the corresponding statistical significance of the observed structure and the errors of the parameters are calculated by doubling the change of likelihood values. The fit with a width fixed to that of the charged structure observed in e + e − → π + π − ψ(3686) [3] yields a mass of M R = (4038.7 ± 6.5) MeV/c 2 (consistent with that of the charged structure M = (4032.1 ± 2.4) MeV/c 2 in Ref. [3]) and a statistical significance of 6.0σ (evaluated by comparing the likelihood values with and without the charmoniumlike state included in the fit). The fit projections on M 2 (π 0 ψ(3686)) and M 2 (π 0 π 0 ) are shown in Fig. 3. Similar to Ref. [3], the fit curves are found to not match the data perfectly. The C.L. of the fit is 19%, estimated by toy-MC studies. An alternative fit with free width yields a mass of M R = (4039.3 ± 6.0) MeV/c 2 , and a width of Γ = (31.9 ± 14.8) MeV, which are consistent with those of the charged structure in Ref. [3] within the statistical uncertainties. Another alternative fit with an additional Z c (3900) 0 included is performed, where the parameters of the Z c (3900) 0 are fixed to the weighted average values M = (3893.6 ± 3.7) MeV/c 2 , Γ = (31.1 ± 7.0) MeV in Refs. [21,23]. The statistical significance of the Z c (3900) 0 is less than 1σ. Similar fits are carried out to the data samples at √ s = 4.258 and 4.358 GeV, respectively, where the parameters of the charmoniumlike state are fixed to those obtained in the data sample at √ s = 4.416 GeV. In the data sample at √ s = 4.258 GeV, a sizable background from e + e − → π + π − ψ(3686) exists. It is included in the fit with the shape fixed to the MC simulation and the magnitude extracted from a fit to the M recoil (π + π − ) spectrum. The statistical significances of the charmoniumlike structure are 3.6σ and 4.5σ for the data samples at √ s = 4.258 and 4.358 GeV, respectively. Alternative fits with additional Z c (3900) 0 states included are performed for the data sample at √ s = 4.258 GeV. Since both Z c (3900) 0 and the structure around 4040 MeV/c 2 are reflected onto each other in the M (π 0 ψ(3686)) spectrum, the statistical significance of Z c (3900) 0 is sensitive to its parameters, and is found to be 1.0σ with the parameters above. The fit procedure has been validated with a set of MC samples. In summary, based on a data sample of e + e − collision data corresponding to 5.2 fb −1 at 16 c.m. energy points between 4.009 and 4.600 GeV collected with the BESIII detector, the Born cross sections for e + e − → π 0 π 0 ψ(3686) at these energy points have been measured for the first time. They are found to be half of those for e + e − → π + π − ψ(3686) [3] within uncertainties, consistent with the expectation from isospin symmetry. The Dalitz plots of π 0 π 0 ψ(3686) are consistent with those in the e + e − → π + π − ψ(3686) [3] at all energy points. Furthermore, a structure is observed in π 0 ψ(3686) with a mass of (4038.7±6.5) MeV/c 2 at √ s = 4.416 GeV, which confirms the structure in the charged mode. No obvious Z c (3900) 0 state is observed in the fit. The new observed structure may provide insight into the properties of the charged structure observed in e + e − → π + π − ψ(3686) as well as the charmoniumlike Z c states observed in analogous decay modes and in charmed meson pairs. However, the fit curve does not match the data perfectly. A future larger statistics sample of data and theoretical input incorporating possible interference effects could lead to a better understanding of the structure.
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center and the supercomputing center of USTC for their strong support. This work is supported in part by National Key