Cross section measurements of $e^+ e^-\to\omega\chi_{c0}$ from $\sqrt{s}=$4.178 to 4.278 GeV

The cross section of the process $e^+e^- \rightarrow \omega \chi_{c0}$ is measured at center-of-mass energies from $\sqrt{s} =$ 4.178 to 4.278 GeV using a data sample of 7 fb$^{-1}$ collected with the BESIII detector operating at the BEPCII storage ring. The dependence of the cross section on $\sqrt{s}$ shows a resonant structure with mass of $(4218.5\pm1.6(\text{stat.})\pm4.0(\text{syst.}))$ MeV/$c^2$ and width of $(28.2\pm3.9(\text{stat.})\pm1.6(\text{syst.}))$ MeV, respectively. This observation confirms and improves upon the result of a previous study. The angular distribution of the $e^+e^- \rightarrow \omega \chi_{c0}$ process is extracted for the first time.

The cross section of the process e + e − → ωχc0 is measured at center-of-mass energies from √ s = 4.178 to 4.278 GeV using a data sample of 7 fb −1 collected with the BESIII detector operating at the BEPCII storage ring. The dependence of the cross section on √ s shows a resonant structure with mass of (4218.5 ± 1.6(stat.) ± 4.0(syst.)) MeV/c 2 and width of (28.2 ± 3.9(stat.) ± 1.6(syst.)) MeV, respectively. This observation confirms and improves upon the result of a previous study. The angular distribution of the e + e − → ωχc0 process is extracted for the first time.
In this paper, we report a study of the e + e − → ωχ c0 reaction based on the most recent e + e − annihilation data collected with the BESIII detector [14] at nine energy points in the range 4.178 √ s 4.278 GeV, with a total integrated luminosity of about 7 fb −1 . The χ c0 state is detected via χ c0 → π + π − /K + K − , and the ω is reconstructed via the ω → π + π − π 0 decay.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [14] located at the Beijing Electron Positron Collider (BEPCII) [15]. 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 the 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the 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. The end cap TOF system was upgraded in 2015 with multigap resistive plate chamber technology, providing a time resolution of 60 ps [16].
Simulated data samples produced with the geant4based [17] 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 models the beam energy spread and initial state radiation (ISR) in the e + e − annihilations using the generator kkmc [18]. For the signal we use a MC sample of the e + e − → ωχ c0 process generated according to the measured angular distribution, which is introduced in section VI. The inclusive MC samples consist of the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [18]. The known decay modes are modelled with evtgen [19] using branching fractions taken from the Particle Data Group (PDG) [20], and the remaining unknown decays from the charmonium states are generated with lundcharm [21]. Final state radiation (FSR) effects from charged final state particles are incorporated via the photos package [22].

III. EVENT SELECTION
For each charged track, the distance of closest approach to the interaction point (IP) is required to be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. The polar angle (θ) of the tracks must be within the fiducial volume of the MDC (| cos θ| < 0.93). Photons are reconstructed from isolated showers in the EMC, which are at least 10 • away from the nearest charged track. The photon energy is required to be at least 25 MeV in the barrel region (| cos θ| < 0.8) or 50 MeV in the end cap region (0.86 < | cos θ| < 0.92). To suppress electronic noise and energy depositions unrelated to the event, the time after the collision at which the photon is recorded in the EMC is required to satisfy 0 ≤ t ≤ 700 ns.
Since the final states of the e + e − → ωχ c0 signal are π 0 π + π − π + π − or π 0 π + π − K + K − , candidate events must have four tracks with zero net charge and at least two photons. The tracks with a momentum larger than 1 GeV/c are identified as π/K from the decay of the χ c0 , whereas lower momentum tracks are considered as pions from ω decays. To reduce the background contributions and to improve the mass resolution, a 5C-kinematic fit is performed to both constrain the total four momentum of the final state partices to the total initial four momentum of the colliding beams and to constrain the invariant mass of the two photons from the decay of the π 0 to its nominal mass [20]. If there is more than one candidate in an event, the one with the smallest χ 2 5C of the kinematic fit is selected. The two track candidates of the decay of the χ c0 are considered to be either a π + π − or a K + K − pair depending on the χ 2 of the 5C-kinematic fit.
, the two tracks are identified as a π + π − pair, otherwise, as a K + K − pair. The χ 2 5C of the candidate events is required to be less than 100.

IV. BORN CROSS SECTION MEASUREMENT
The correlation between the π + π − π 0 invariant mass, M (π + π − π 0 ), and the π + π − /K + K − mass, M (π + π − /K + K − ), is shown in the top panel in Fig. 1 for data taken at √ s = 4.219 GeV. A high density area can be observed that originates from the e + e − → ωχ c0 channel. The mass range [0.75, 0.81] GeV/c 2 in M (π + π − π 0 ) is defined as the ω signal region and is in-dicated by horizontal dashed lines. A sideband in the range [0.60, 0.72] GeV/c 2 is used to study the nonresonant background. The χ c0 signal region is indicated by the vertical dashed lines and is defined as [3.38, 3.45] GeV/c 2 . The bottom panel of Fig. 1 shows the distribution of M (π + π − /K + K − ) for data in the ω signal region. The shaded (green) histogram corresponds to normalized events in the ω sideband region.  An unbinned maximum likelihood fit is performed to obtain the signal yields. In the fit, we use the MCdetermined shape to describe the χ c0 signal. The background is described with a generalized ARGUS function [23] where m 0 is fixed to ( √ s − 0.75 GeV), with 0.75 GeV being the lower limit of M (π + π − π 0 ), and p, k are free parameters. The red solid curve in the bottom panel of Fig. 1 shows the fit result. The data taken at the other center-of-mass energies are analyzed using the same method and the fit results are summarized in Table I. for the e + e − → ωχc0 reaction at the different center-ofmass energies √ s, together with integrated luminosities L, number of signal events N sig , radiative correction factor 1 + δ(s), vacuum polarization factor 1 |1−Π| 2 , and efficiency ǫ. The first uncertainties are statistical, and the second systematic. The Born cross section is calculated with where N sig is the number of signal events, L is the integrated luminosity obtained using the same method in Ref. [24], 1 + δ(s) is the radiative correction factor obtained from a Quantum Electrodynamics (QED) calculation [18,25] using the obtained preliminary cross section as input and iterating it until the results converge, 1 |1−Π| 2 is the correction factor for vacuum polarization [26], B is the product of branching fractions B(χ c0 → π + π − /K + K − ) × B(ω → π + π − π 0 ) × B(π 0 → γγ), and ǫ is the event selection efficiency. The Born cross sections (or upper limits at 90% C.L.) at each energy point for e + e − → ωχ c0 are listed in Table I.
The luminosity is measured with a precision of about 1.0% using the well-known Bhabha scattering process [24]. The uncertainty in the tracking efficiency is obtained as 1.0% per track using the process e + e − → π + π − K + K − [2]. The uncertainty in photon reconstruction is 1.0% per photon, obtained by studying the J/ψ → ρ 0 π 0 decay [27].
The systematic uncertainty due to the kinematic fit is estimated by correcting the helix parameters of charged tracks according to the method described in Ref. [28]. The difference between detection efficiencies obtained from MC samples with and without correction is taken as the uncertainty.
The line-shape of the e + e − → ωχ c0 cross section will affect the radiative correction factor and the efficiency. In the nominal results, we use a piecewise Gaussian function as the line-shape to describe the cross section. The shape is used as input and is iterated until the results converge. To estimate the uncertainty from the radiative correction, we change the line-shape to the published Breit-Wigner (BW) function of the Y (4220) [2]. The difference between the results is taken as a systematic uncertainty.
The uncertainty from the fit range is obtained by varying the limits of the fit range by ±0.01 GeV/c 2 . We take the largest difference of the corresponding cross section measurement with respect to the nominal one as the systematic uncertainty. We use the MC-determined shape convolved with a Gaussian function to fit the data as input to get the uncertainty of the signal shape. The difference in the results with respect to the nominal one is taken as the systematic uncertainty. To estimate the systematic uncertainty caused by the background shape, we vary m 0 by ±0.01 GeV/c 2 in the ARGUS function, and take the largest difference in the results as the uncertainty.
The measured angular distribution is used as a model to generate signal events in the MC simulations. The detection efficiency of the e + e − → ωχ c0 reaction will depend upon its angular distribution. We obtained an angular distribution parameter, defined in section VI, of α = −0.30 ± 0.18(stat.) ± 0.05(sys.). The systematic uncertainty of the efficiency due to uncertainties in the angular distribution is estimated by varying the α value by one standard deviation, the total uncertainty on α.
The uncertainty in the product of the branching fractions B(χ c0 → π + π − /K + K − ) × B(ω → π + π − π 0 ) × B(π 0 → γγ) is taken from the uncertainties quoted by the PDG [20]. Table II summarizes all the systematic uncertainties related to the cross section measurements of the e + e − → ωχ c0 process for each center-of-mass energy. The overall systematic uncertainties are obtained by adding all the sources of systematic uncertainties in quadrature assuming they are uncorrelated.  Figure 2 shows the dressed cross sections (σ = σ B |1−Π| 2 ) for the e + e − → ωχ c0 reaction as a function of centerof-mass energy. The black square points are taken from Refs. [1,2], and the blue circular points are from this work. We observe an enhancement in the cross section around 4.22 GeV. By assuming that the ωχ c0 signals all come from a single resonance, which we label as the Y (4220), with mass M and width Γ, we fit the cross section data with the following formula convolved with a Gaussian function for the energy spread: where Φ( √ s) is the two-body phase space factor and Γ ee is the electronic width. The fit to all the data in Fig. 2 gives Γ ee B(ωχ c0 ) = (2.5 ± 0.2) eV, M = (4218.5 ± 1.6) MeV/c 2 , Γ = (28.2 ± 3.9) MeV, where the uncertainties are statistical only. In the fit, the cross sections' statistical uncertainties are used only. The goodness of fit is χ 2 /ndf = 29/19, where ndf is the number of degrees of freedom.
The systematic uncertainties on the resonant parameters mainly arise from uncertainties in the absolute beam energy, the parametrization of the BW function, and the cross section measurement. Because the energy spread effect has been considered in the fit, we ignore the systematic uncertainty from energy spread.
Since the uncertainty of the beam energy is about 0.8 MeV, which is obtained using the same method in Ref. [29], the uncertainty of the resonant parameters caused by the beam energy is estimated by varying √ s within 0.8 MeV.
The cross section has been fitted with a BW function having the energy-dependent width Γ = Γ 0 Φ( √ sult and the nominal result is taken as the uncertainty from the parametrization of the BW function. The systematic uncertainty of the cross section measurement will affect the resonant parameters in the fit and can be divided into two parts. One part comes from the uncorrelated uncertainty among the different centerof-mass energies, and the other part is a common uncertainty. The first part has been considered by including the systematic uncertainty of the cross section in the fit. The difference between the parameters obtained in this fit to those from the nominal fit is taken as the uncertainty. We vary the cross section within the systematic uncertainty coherently for the second part and take the difference between this fit result and the nominal result as the uncertainty. We add the two parts in quadrature assuming they are uncorrelated.
Table III summarizes all the systematic uncertainties of the resonant parameters. The total systematic uncertainty is obtained by summing all the sources of sys-tematic uncertainties in quadrature by assuming they are uncorrelated.

VI. ANGULAR DISTRIBUTION MEASUREMENT
Both S and D-wave contributions are possible in the process Y (4220) → ωχ c0 . A measurement of their strengths can be helpful to extract information about the underlying dynamics of the decay process. We therefore performed an angular analysis of the relatively highstatistics data samples taken at √ s = 4.219, 4.226, and 4.236 GeV (selection of √ s = 4.226 GeV data was reported in Ref. [1]). The helicity angle, θ ω , defined by the scattering angle of the ω with respect to the electron beam in the e + e − center-of-mass frame was reconstructed for each event. Figure 3 shows the bin-by-bin efficiency-corrected events as a function of cosθ ω for the three center-of-mass energies. We performed a simultaneous fit using the function 1+αcos 2 θ ω , assuming α is common to the three energy points. The red line in Fig. 3 shows the best fit result with α = −0.30 ± 0.18 ± 0.05, where the first uncertainty is statistical and the second systematic. The goodness of the fit is χ 2 /ndf= 31/26. The fit result indicates evidence for a combination of S and D−wave contributions in the Y (4220) → ωχ c0 process, although the statistical significance of this conclusion is only 2σ compared with a pure S−wave contribution of α = 0.
The systematic uncertainty of α has been estimated by varying the fit range (0.02), the signal (0.01) and background shapes (0.03), and the radiative correction factor (0.03). The uncertainties are indicated in brackets and are determined with the same method described earlier for the cross section measurements. In addition, we estimate an additional source of systematic uncertainty by varying the number of bins. For this, we change the number of bins from 10 to 8, and repeat the process. The difference in α is found to be 0.01. The overall systematic uncertainty (0.05) is obtained by summing all the items of systematic uncertainties in quadrature by assuming they are uncorrelated. shows the summed result of the three center-of-mass energies.

VII. SUMMARY
The process e + e − → ωχ c0 has been studied using 9 data samples collected at center-of-mass energies from √ s = 4.178 to 4.278 GeV. The √ s-dependence of the cross section has been measured and the results are listed in Table I and are shown in Fig. 2. A clear enhancement is seen around √ s = 4.22 GeV which confirms, and statistically improves upon, an earlier observation [1]. By fitting the e + e − → ωχ c0 cross section with a single resonance, the mass and width for the structure are determined to be M = (4218.5±1.6(stat.)±4.0(sys.)) MeV/c 2 and Γ = (28.2±3.9(stat.)±1.6(sys.)) MeV. The obtained resonance parameters are not compatible with the vector charmonium state ψ(4160), ruling out its possible contribution to the structure [6]. Moreover, we studied the angular distribution of the process Y (4220) → ωχ c0 . We measured α = −0.30±0.18±0.05, which indicates a combination of S and D−wave contributions in the decay. Figure 4 shows the measured mass and width of the Y (4220) from the different processes. The masses are consistent with each other, while the widths are not. The widths from the processes e + e − → π + π − h c , π + π − ψ(3686), and π + D 0 D * − + c.c. are larger than those from the processes e + e − → ωχ c0 and π + π − J/ψ. From these inconsistencies in the width, we cannot draw a conclusion on whether the structure observed in these processes is the same state or whether the inconsistencies are caused by the BW parameterization. Further experimental studies with higher statistics are needed to draw a more reliable conclusion on the nature of this structure. Mass and width of the Y (4220) obtained from the processes e + e − → ωχc0, π + π − hc, π + π − J/ψ, π + π − ψ(3686) and π + D 0 D * − + c.c.