Study of the process $e^+e^-\to\pi^0\pi^0 J/\psi$ and neutral charmonium-like state $Z_c(3900)^0$

Cross sections of the process $e^+e^-\to\pi^0\pi^0 J/\psi$ at center-of-mass energies between 3.808 and 4.600 GeV are measured with high precision by using 12.4 $fb^{-1}$ of data samples collected with the BESIII detector operating at the BEPCII collider facility. A fit to the measured energy-dependent cross sections confirms the existence of the charmonium-like state $Y(4220)$. The mass and width of the $Y(4220)$ are determined to be $(4220.4\pm2.4\pm2.3)$ ${\rm MeV}/c^{2}$ and $(46.2\pm4.7\pm2.1)$ MeV, respectively, where the first uncertainties are statistical and the second systematic. The mass and width are consistent with those measured in the process $e^+e^-\to\pi^+\pi^- J/\psi$. The neutral charmonium-like state $Z_c(3900)^0$ is observed prominently in the $\pi^0 J/\psi$ invariant-mass spectrum, and, for the first time, an amplitude analysis is performed to study its properties. The spin-parity of $Z_c(3900)^0$ is determined to be $J^{P}=1^{+}$, and the pole position is $(3893.1\pm2.2\pm3.0)-i(22.2\pm2.6\pm7.0)$ ${\rm MeV}/c^2$, which is consistent with previous studies of electrically charged $Z_c(3900)^\pm$. In addition, cross sections of $e^+e^-\to\pi^0 Z_c(3900)^0\to \pi^0\pi^0 J/\psi$ are extracted, and the corresponding lineshape is found to agree with that of the $Y(4220)$.

It is particularly interesting to understand the relationship between the Y and Z c production through the same decay processes [3,7,8,14,15], since such information could shed light on their nature. More precise measurements of their resonant parameters, production cross sections, and decay modes, as well as searches for new states, are essential. Studying the neutral decay process e + e − → π 0 π 0 J/ψ is a natural way to study the properties of the Y (4220) and to search for the chargeneutral isospin partner of the Z c (3900) ± .
The neutral Z c (3900) 0 was observed in the processes e + e − → π 0 π 0 J/ψ and e + e − → π 0 (DD * ) 0 by CLEO-c and BESIII [17][18][19], demonstrating that the Z c (3900) ± is an isovector. However, the previous analysis was unable to extract the resonant parameters of the Y states and to determine the properties of the Z c (3900) 0 , due to limited statistics and lack of data in some key center-of-mass (c.m.) energy points. In this Letter, we present an updated analysis of e + e − → π 0 π 0 J/ψ using 12.4 fb −1 [20] of BESIII data taken at c.m. energies from 3.808 to 4.600 GeV. The production cross sections are measured with high precision, and a PWA is performed to study properties of the Z c (3900) 0 .
Candidates for charged tracks and photons are selected using the same criteria as described in Ref. [18]. Electrons (muons) identification and π 0 reconstruction are proceeded with the same strategies as reported in Refs. [24,25], respectively. Less than three π 0 π 0 combinations for each event are required to suppress combinatorial backgrounds. A kinematic fit under the hypothesis e + e − → π 0 π 0 + − is applied to enforce energymomentum conservation and to constrain the two π 0 masses. The combination with the smallest χ 2 6C and with χ 2 6C < 75 is chosen. Detailed MC studies indicate that the dominant backgrounds are those from e + e − → qq, e + e − → ηJ/ψ (η → π 0 π 0 π 0 ) and γψ(2S) (ψ(2S) → π 0 π 0 J/ψ). The first one has a uniform distribution in the + − invariantmass (M + − ) spectrum, while the other two produce a broad peak around the J/ψ position. To determine the signal yield of e + e − → π 0 π 0 J/ψ, an unbinned maximumlikelihood fit to the M + − spectrum in the range 2.90 < M + − < 3.30 GeV/c 2 is performed. The signal and peaking backgrounds are described using MC-simulated shapes convolved with a common Gaussian function representing the resolution difference between data and the MC simulation, while the other possible backgrounds are described with a first-order polynomial function. The peaking background contributions are normalized to the integrated luminosity.
The Born cross section is determined by where N obs is the signal yield obtained as described above, L int is the integrated luminosity, ε is the selection efficiency estimated with signal MC samples modeled ac-cording to the results of the PWA for those points that fall within √ s = 4.178 ∼ 4.416 GeV and by PHSP for other energy points with low statistics, (1 + δ r ) is the radiative correction factor, (1 + δ v ) is the vacuum polarization factor derived from QED calculations [26] and B is the world averaged value [27] of the product branching fractions (BFs) of the intermediate states. The resultant cross sections are shown in the top subplot of Fig. 1, and the various numbers used in the calculation are summarized in Table I of the supplementary material [20]. To examine the isospin symmetry, the cross-section ratios of the processes e + e − → π 0 π 0 J/ψ to e + e − → π + π − J/ψ [4] at different c.m. energies are extracted, and are fitted by a constant function, as shown in the bottom subplot of Fig. 1. The fit yields an average ratio of 0.48 ± 0.02, consistent with the prediction (0.5) based on isospin symmetry. (GeV) s  4220) and Y (4320), respectively. Bottom: Cross-section ratio of e + e − → π 0 π 0 J/ψ to e + e − → π + π − J/ψ, where the blue dashed line corresponds to the average.
A χ 2 fit is performed to the Born cross section to study resonant structure in the range √ s = 3.808 ∼ 4.600 GeV. The probability density function (PDF) is defined as represents the non-resonant component [28], M thd = 2m π 0 + m J/ψ , Φ( √ s) is the PHSP factor of the three-body decay R i → π 0 π 0 J/ψ [27], f i stands for the amplitude of structure R i (R 1 is Y (4220) and R 2 is Y (4320)), and φ i is its phase relative to the continuum. f i is defined as where, M i , Γ i tot , Γ i ee and B i π 0 π 0 J/ψ are the mass, full width, partial width coupling to e + e − , and the decay BF of R i → π 0 π 0 J/ψ, respectively. The mass and width of the Y (4320) are fixed to those reported in Ref. [4] and other parameters are allowed to vary. The fit yields four solutions with equal fit qualities (the goodness of fit is χ 2 /ndf = 19.3/19, where ndf is the number of degrees of freedom), and (4220.4 ± 2.4) MeV/c 2 and (46.2 ± 4.7) MeV for the mass and width of the Y (4220), respectively, which agree with the ones obtained in e + e − → π + π − J/ψ [4]. The fit curves for one of the solutions are shown in Fig. 1. A summary of the fit results can be found in Table II of the supplementary material [20]. The statistical significance of Y (4320) is estimated to be 4.2σ by the changes in the likelihood values obtained from including and excluding the Y (4320).
The systematic uncertainties for the cross section measurement include those associated with the luminosity, detection efficiency, radiative correction, fitting procedure, and the BFs of intermediate states. The uncertainties of the detection efficiency include those of tracking, photon detection, e/µ identification, π 0 reconstruction, kinematic fit, and MC modeling. The uncertainties for the Y (4220) resonant parameters are from the measured cross section, the c.m. energy measurement and spread, and the fit procedure. All above uncertainties are estimated, and a detailed description can be found in the supplementary material [20]. Assuming all sources are independent, the total uncertainties are obtained by adding all the values of the individual contributions in quadrature, resulting in 6.1% for the Born cross section, 2.3 MeV/c 2 and 2.1 MeV for the mass and width of Y (4220), respectively.
For the PWA, candidate signal events are required to be within the J/ψ signal region (3.06 < M + − < 3.13 GeV/c 2 ). Events within the J/ψ-sideband region ((2.93, 3.00) ∪ (3.20, 3.27)) GeV/c 2 are used to evaluate the effect of non-J/ψ background. The J/ψ-peaking background is found to be negligible. To improve the purity of samples in the J/ψ → µ + µ − channel, at least one muon track penetrating more than five muon chamber (MUC) layers is required. The invariant-mass distributions of π 0 J/ψ (M π 0 J/ψ ) and π 0 π 0 (M π 0 π 0 ) for the four data samples of √ s = 4.226, 4.236, 4.244 and 4.258 GeV, which have the largest statistics among all data samples, are shown in Fig. 2, where the structures Z c (3900) 0 and f 0 (980) are clearly visible. To decompose the intermediate processes, a simultaneous PWA is performed on the four data samples with a sum of 2188 ± 47 candidate events, of which 57 ± 8 are background events estimated from the J/ψ-sideband regions.
The isobar model is implemented in the PWA and includes the cascading decay chains e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ and e + e − → RJ/ψ → π 0 π 0 J/ψ, where R represents σ, f 0 (980), f 0 (1370), etc.. The quasi two-body decay amplitudes in the sequential decays are constructed using the helicity formalism [29]. The amplitude of J/ψ → + − is considered, and an additional term is introduced to correct for the difference in the azimuthal angle between the lepton helici- Points with errors are data, red solid curves are the total fit results, the blue dashed (magenta long-dashed) curves represent Zc(3900) 0 (π 0 π 0 -S wave) components, and green shaded histograms represent the estimated backgrounds. Each event appears twice in the M π 0 J/ψ distributions. The χ 2 /ndf is calculated by merging those bins with less than 10 events. ties in two different decay chains [30]. The intermediate states are parameterized with the relativistic Breit-Wigner (BW) functions, except for the f 0 (980), which is described by a Flatté formula with fixed parameters taken from Ref. [31]. The width of the wide σ resonance is parameterized as 1 − 4m 2 π s Γ [32]. The resonant parameters of the f 0 (1370) and other well-known mesons are taken from the world averaged values [27], while those of the Z c (3900) 0 are left free in the fit. The relative magnitudes and phases of the individual intermediate processes are determined by performing an unbinned maximum-likelihood fit for the four data samples simultaneously. The properties of the intermediate states, i.e., their masses, widths, and coupling strength to the different final states are shared parameters between the different data samples, while the production magnitudes and the relative phases are independent. The background contribution to the overall likelihood value is estimated from the events in the J/ψ-sideband region and subtracted. All the possible known intermediate states are considered in the fit and only those with statistical significances larger than 5σ are kept in the nominal results. The statistical significance is calculated by the differences in the likelihood values and the number of degrees of freedom between the two scenarios with or without including the intermediate process. Different Z c (3900) 0 spin-parity hypotheses are tested in the simultaneous fit.
With the above strategy, the nominal fit includes the intermediate processes e + e − → σJ/ψ, f 0 (980)J/ψ, f 0 (1370)J/ψ and π 0 Z c (3900) 0 . As shown in Fig. 2, the π 0 π 0 S-wave contribution dominates. The spin-party of the Z c (3900) 0 is determined to be J P = 1 + with a statistical significance of more than 9σ over alternative J P hypotheses (J P = 0 + , 1 − , 2 + , 2 − ). The fits yield a mass (3893.0 ± 2.3) MeV/c 2 and a width (44.2 ± 5.4) MeV for the Z c (3900) 0 , where the uncertainties are statistical only. Since the Z c (3900) 0 is also observed in e + e − → π 0 D * D [19] and its mass is close to the mass threshold of D * D , we also perform an alternative fit by parameterizing the Z c (3900) 0 with a Flatté formula as those for e + e − → π + π − J/ψ [16]. The fit results are not sensitive enough to determine the coupling constants g π 0 J/ψ and g D * D . However, if the ratio g D * D /g π 0 J/ψ is fixed to the value reported in Ref. [16], the fit results in a comparable description of the data as the nominal case.
Based on the above procedure, the Born cross sections for the process e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ are measured using σ Born c is the fraction of the Z c (3900) 0 component in the e + e − → π 0 π 0 J/ψ process, extracted from the PWA. To extract the energy-dependent lineshape of the cross section for e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ around the Y (4220), we also perform the PWA for other data samples with c.m. energy around the Y (4220), individually, where the properties of the Z c (3900) 0 and its coupling strength to π 0 J/ψ are fixed to those obtained from the simultaneous fit. The resulting cross sections, as shown in Fig. 3 and in Table V of the supplementary material [20], show a clear structure around 4220 MeV. We perform a χ 2 fit to the Z c (3900) 0 cross sections using a coherent sum of a relativistic BW function and σ N Y ( √ s), as in formula 2, but with a threshold M thd = M π 0 + M Zc(3900) 0 for the non-resonant component and a PHSP factor of two-body decay R → π 0 Z c (3900) 0 . The fit curve is shown in Fig. 3, with a goodness of fit of χ 2 /ndf = 8.5/5. The mass (4231.9 ± 5.3) MeV/c 2 and width (41.2 ± 16.0) MeV of the resonant structure are consistent with those of the Y (4220) presented previously. We also tested the scenarios including the Y (4320) with fixed resonant parameters and/or phase, but none of those improve the fit quality. Due to the lack of data around 4.3 GeV, we cannot rule out a contribution from the Y (4320).
The systematic uncertainties for the Z c (3900) 0 resonant parameters and the corresponding cross sections include those associated with the amplitude modeling and background treatment in the PWA, as well as the detection efficiency difference between data and the MC simulation. The uncertainties associated with the amplitude modeling in PWA stem from the parameterizations of intermediate states (σ, f 0 (980), f 0 (1370) and Z c (3900) 0 ), the barrier radius in the barrier factor, and possible missing components. We perform alternative PWAs for the scenarios with different parameterizations of the intermediate states, different barrier radius, and an additional e + e − → f 2 (1270)J/ψ amplitude individually, and the resultant (largest) changes with respect to the nominal values are regarded as the uncertainties. The uncertainty related to the background treatment is studied by varying the M + − -sideband region. The uncertainties related with the detection efficiency include those for tracking of leptons and photons, π 0 reconstruction, kinematic fit and ISR correction factor, as well as the uncertainty by requiring additional MUC hits. The detection resolution, about 8.8 MeV/c 2 in the M π 0 J/ψ distribution, is not considered in the nominal PWA. Its effect is estimated with 300 sets of pseudo-experiments. All the above uncertainties are summarized in Table VI of the supplemental material [20] including a more detailed description.
The uncertainties for the resonant parameters in the fit to σ(e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ) include the same sources as discussed in the fit to σ(e + e − → π 0 π 0 J/ψ) and the same approaches are used, as summarized in Table VII of the supplemental material [20]. In this case, the uncertainty of the fit procedure is estimated by replacing the nominal PDF model with one BW function only and by changing the constant full width of BW to a phase-space dependent width, Γ Φ( In summary, we measured Born cross sections of e + e − → π 0 π 0 J/ψ for c.m. energies between 3.808 and 4.600 GeV with data samples collected by the BESIII experiment. The measured cross sections are fitted by including two resonant structures, Y (4220) and Y (4320), with the resonant parameters of the Y (4320) fixed to the values taken from Ref. [4]. The mass and width of the Y (4220) are measured to be (4220.4 ± 2.4 ± 2.3) MeV/c 2 and (46.2 ± 4.7 ± 2.1) MeV, respectively, where the first uncertainties are statistical, and the second are systematic (the same as following). These measurements agree with those reported in Ref. [4], and confirm the existence of the Y (4220). The average ratio of the cross section e + e − → π 0 π 0 J/ψ to that of e + e − → π + π − J/ψ [4] is 0.48 ± 0.02, which is consistent with isospin symmetry.
The Z c (3900) 0 signal is clearly observed in the M π 0 J/ψ distribution, and a PWA is performed to study its properties. The spin-parity of the Z c (3900) 0 is determined to be J P = 1 + , and the measured mass (3893.0 ± 2.3 ± 19.9) MeV/c 2 and width (44.2 ± 5.4 ± 9.1) MeV correspond to a pole position (3893.1±2.2±3.0)−i(22.2±2.6± 7.0) MeV/c 2 . These values are consistent with those of the charged Z c (3900) ± observed in e + e − → π + π − J/ψ. The Born cross sections of e + e − → π 0 Z c (3900) 0 → π 0 π 0 J/ψ are also measured and fitted for c.m. energies between 4.178 and 4.278 GeV. The fit yields a structure with a mass of (4231.9 ± 5.3 ± 4.9) MeV/c 2 and a width of (41.2 ± 16.0 ± 16.4) MeV, compatible with the Y (4220). The results indicate a strong correlation between the production of the Y (4220) and the Z c (3900) 0 . Due to the lack of data around 4.3 GeV, the existence of the Y (4320) in the Z c (3900) 0 production cannot be ruled out.