Measurement of $e^+e^-\to\gamma \chi_{c0,c1,c2}$ cross sections at center-of-mass energies between 3.77 and 4.60 GeV

The $e^+e^-\to\gamma\chi_{cJ}$ ($J=0,1,2$) processes are studied at center-of-mass energies ranging from 3.773 to 4.600 GeV, using a total integrated luminosity of 19.3 fb$^{-1}$ $e^+e^-$ annihilation data accumulated with the BESIII detector at BEPCII. We observe for the first time $e^+e^-\to\gamma \chi_{c1,c2}$ signals at $\sqrt{s}=$ 4.180 GeV with statistical significances of 7.6$\sigma$ and 6.0$\sigma$, respectively. The production cross section of $e^+e^-\to\gamma \chi_{c1,c2}$ at each center-of-mass energy is also measured. We find that the line shape of the $e^+e^-\to\gamma \chi_{c1}$ cross section can be described with conventional charmonium states $\psi(3686)$, $\psi(3770)$, $\psi(4040)$, $\psi(4160)$. Compared with this, for the $e^+e^-\to\gamma \chi_{c2}$ channel, one more additional resonance is added to describe the cross section line shape. Its mass and width are measured to be $M=4371.7\pm7.5\pm1.8$ MeV/$c^2$ and $\Gamma^{tot}=51.1\pm17.6\pm1.9$ MeV, where the first uncertainties are statistical and the second systematic. The significance of this resonance is estimated to be 5.8$\sigma$, and its parameters agree with the $Y(4360)$ resonance previously reported in $e^+e^-\to\pi^+\pi^-\psi(3686)$, and the $Y(4390)$ in $e^+e^-\to\pi^+\pi^-h_c$ within uncertainties. No significant signal for the $e^+e^-\to\gamma \chi_{c0}$ process is observed, and the upper limits of Born cross sections $\sigma_{B}(e^+e^-\to\gamma \chi_{c0})$ at 90\% confidence level are reported.

The e + e − → γχcJ (J = 0, 1, 2) processes are studied at center-of-mass energies ranging from 3.773 to 4.600 GeV, using a total integrated luminosity of 19.3 fb −1 e + e − annihilation data accumulated with the BESIII detector at BEPCII. We observe for the first time e + e − → γχc1,c2 signals at √ s = 4.180 GeV with statistical significances of 7.6σ and 6.0σ, respectively. The production cross section of e + e − → γχc1,c2 at each center-of-mass energy is also measured. We find that the line shape of the e + e − → γχc1 cross section can be described with conventional charmonium states ψ(3686), ψ(3770), ψ(4040), ψ(4160). Compared with this, for the e + e − → γχc2 channel, one more additional resonance is added to describe the cross section line shape. Its mass and width are measured to be M = 4371.7 ± 7.5 ± 1.8 MeV/c 2 and Γ tot = 51.1 ± 17.6 ± 1.9 MeV, where the first uncertainties are statistical and the second systematic. The significance of this resonance is estimated to be 5.8σ, and its parameters agree with the Y (4360) resonance previously reported in e + e − → π + π − ψ(3686), and the Y (4390) in e + e − → π + π − hc within uncertainties. No significant signal for the e + e − → γχc0 process is observed, and the upper limits of Born cross sections σB(e + e − → γχc0) at 90% confidence level are reported.

I. INTRODUCTION
In the past decades, many charmoniumlike states were observed experimentally, such as the X(3872), Y (4260), and Z c (3900) [1]. Among them, the vector Y -states should have quantum numbers J P C = 1 −− , as they are produced in e + e − annihilation process. Considering the Y (4260) [2][3][4][5], Y (4360), and Y (4660) states [6][7][8][9], together with the conventional charmonium states ψ(4040), ψ(4160), and ψ(4415), there are at least six vector states between 4.0 and 4.7 GeV. However, the potential model only predicts five vector charmonium states in this mass region [10]. In addition, unlike the known 1 −− conventional charmonium states that decay predominantly into open-charm final states [D ( * ) D ( * ) ], the Y -states show strong coupling to hidden-charm final states. These unusual behaviors indicate that the Y -states might be nonconventional quarkonium states. To better understand the nature of these states and also one gets better insights in the relevant degrees of freedom that play a role in these systems that are governed by the strong interaction, it is important to further investigate these states experimentally.
Experimentally, the e + e − → γχ cJ (J = 1, 2) processes above 4 GeV have been studied before by BESIII [15], CLEO [16], and Belle experiments [17]. Due to the limited statistics, no obvious signal has been observed between 4 -5 GeV. The BESIII has collected the world's largest dataset from 4.0 to 4.6 GeV, and it is thus highly motivated to search for these decay modes.
In this paper, we report the study of the e + e − → γχ cJ (J = 0, 1, 2) processes at e + e − center-of-mass (c.m.) energies between √ s = 4.008 -4.6 GeV, using data samples corresponding to an integrated luminosity of 16.0 fb −1 accumulated with the BESIII detector at the BEPCII collider. To better estimate the contributions from ψ(3686) and ψ(3770), the datasets with integrated luminosity of 3.3 fb −1 between √ s of 3.773 and 4.008 GeV for e + e − → γχ c1,c2 channels are also analyzed. The datasets together with the corresponding c.m. energies are summarized in Table A in the Appendix. Compared with the previous BESIII measurement [15], the new dataset covers an extended c.m. energies with about one order of magnitude higher luminosity, and also both J/ψ → e + e − /µ + µ − events (only µ + µ − used in previous work) are studied. The integrated luminosities are measured with Bhabha events (e + e − → (γ)e + e − ) with an uncertainty of 1% [18]. The c.m. energy of each dataset is measured using dimuon events (e + e − → (γ)µ + µ − ), with an uncertainty of ± 0.8 MeV [19].

II. BESIII DETECTOR AND MC SIMULATION
The BESIII detector [20] records symmetric e + e − collisions provided by the BEPCII storage ring [21], with a designed peak luminosity of 1 × 10 33 cm −2 s −1 at c.m. energy of 3.77 GeV. BESIII has collected large data samples between 2.0 and 4.6 GeV [22]. The cylindrical core of the BESIII detector covers 93% of the full solid angle and 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 identification modules interleaved with steel. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for 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 [20]. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [23].
Simulated Monte Carlo (MC) samples produced with geant4-based [24] software, which includes the geometrical description of the BESIII detector and the detector response, are used to determine the detection efficiency, and to estimate physical background. The signal MC e + e − → γχ c0,c1,c2 events are generated assuming a pure E1 transition. The simulation models the beam energy spread and initial-state-radiation (ISR) in e + e − annihilation using the generator kkmc [25]. The maximum ISR photon energy is set to the energy corresponding to the γχ c0,c1,c2 production threshold. The final-stateradiation (FSR) from charged final state particles is modelled with photos [26]. Possible background contributions are investigated with the inclusive MC samples, which consist of open-charm processes, the ISR production of lower mass vector charmonium(-like) states, and the continuum processes. The known decay modes of charmed hadrons are modelled with evtgen [27], with known branching fractions taken from the Particle Data Group (PDG) [12], and the remaining unknown decays with lundcharm [28].

A. Event selection
The final state particles for e + e − → γχ c1,c2 are γγℓ + ℓ − , where the χ c1,c2 are reconstructed with γJ/ψ, and the J/ψ is reconstructed with ℓ + ℓ − (ℓ = e or µ). Events with two charged tracks with zero net charge and at least two photons are selected. Each charged track is required to originate from the interaction point, within ±1 cm in the plane perpendicular to the beams and 10 cm along the beam direction. The |cos θ| of each charged track is required to be less than 0.93, where θ is the polar angle of each track. Photons are required to have a deposited energy larger than 25 MeV in the barrel EMC region (|cos θ| < 0.8) and larger than 50 MeV in the end cap region (0.86 < |cos θ| < 0.92). The EMC time for a photon is required to be within 700 ns of the event start time to suppress the electronic noise and energy deposition unrelated to the physical events. Each charged track should have a momentum larger than 1 GeV/c. For leptons, we use the energy deposited in the EMC to separate electrons from muons. Charged tracks with the energy deposited in the EMC larger than 1 GeV are identified as electrons, and charged tracks with the energy deposited less than 0.4 GeV are identified as muons. For photons, the two most energetic photons are regarded as the candidates for signal events. Through the paper, we denote the photon with higher energy as γ H , and the other as γ L .
To improve the mass resolution and to suppress backgrounds, a four-constraint (4C) kinematic fit is performed under the γγℓ + ℓ − hypothesis, which constrains the total four momentum of the final measured particles to the initial four-momentum of the colliding beams. The χ 2 of the kinematic fit is required to be less than 40.
To suppress the radiative Bhabha events (e + e − → γe + e − ) in the J/ψ → e + e − mode, the cosine of the opening angle between the electron and the nearest photon (cos θ eγ ) is required to be less than 0.86. Since the photon from radiative Bhabha process is always close to the beam direction, the cosine of the polar angle of the selected photons are required to satisfy |cos θ γ L/H | < 0.8. In both e + e − and µ + µ − modes, the background from e + e − → ηJ/ψ with η → γγ is rejected by requiring MC simulations show that the background from e + e − → γ ISR ψ(3686) with ψ(3686) → γχ c1,c2 can be ignored for most of the energies expect for data at √ s = 3.773 GeV. These background events are simulated and subtracted from the signal yield at √ s = 3.773 GeV. A fit to the lepton pair invariant mass gives a resolution of 10.8 MeV/c 2 and 10.5 MeV/c 2 for J/ψ → e + e − and J/ψ → µ + µ − events, respectively. The J/ψ mass window is defined as 3.08 < M (ℓ + ℓ − ) < 3.12 GeV/c 2 . While the sidebands of the J/ψ are defined by 3.00 < M (ℓ + ℓ − ) < 3.06 GeV/c 2 and 3.14 < M (ℓ + ℓ − ) < 3.20 GeV/c 2 , which is three times as wide as the J/ψ signal region.

B. Cross section
According to kinematic, we find that the photon from e + e − → γχ c1,c2 has lower energy than the one from χ c1,c2 → γJ/ψ for data with √ s < 4.009 GeV. On the contrary, the former has higher energy than the latter for data with √ s > 4.009 GeV. To obtain the number of signal events, we make use of both fitting and counting methods. For each data sample with L int > 400 pb −1 (L int is integrated luminosity), a fit with J/ψ → e + e − or µ + µ − events is performed to the invariant mass distribution of γ H J/ψ ( √ s < 4.009 GeV) or these two photons cannot be distinguished by energy. A 2-dimensional fit to the distribution of M (γ H J/ψ) versus M (γ L J/ψ) is used to extract the number of signal events. In these fits, the signal probability density functions (PDFs) are described with MC-simulated shapes, and the background PDFs are constrained to J/ψ sideband events. For the low-statistics data samples with L int < 200 pb −1 and √ s > 4.009 GeV, we obtain the signal yield by counting the number of events in the χ c1,c2 signal region and by subtracting the number of normalized background events in the χ c1,c2 mass sideband region. The χ c1 and χ c2 signal regions are defined as 3.49 < M (γ L J/ψ) < 3.53 GeV/c 2 and 3.54 < M (γ L J/ψ) < 3.58 GeV/c 2 , which include more than 94% of the signal events. The sidebands of χ c1,c2 are defined as 3.42 < M (γ L J/ψ) < 3.46 GeV/c 2 and 3.6 < M (γ L J/ψ) < 3.64 GeV/c 2 . Taking √ s = 4.178 GeV as an example, the invariant mass distribution of M (γ L J/ψ) as well as the fit results for the surviving events are shown in Fig. 1. Clear χ c1,c2 signals are observed. The statistical significance of χ c1,c2 signals are calculated by comparing the log-likelihoods with and without the signal components in the fit, and taking the change of number of degrees of freedom into account. The statistical significances are estimated to be 7.6σ for the χ c1 signal and 6.0σ for the χ c2 signal. This is the first observation of the e + e − → γχ c1,c2 processes between 4 -5 GeV. The invariant mass distributions of The production cross section of e + e − → γχ c1,c2 at each e + e − c.m. energy is calculated as where N signal is the number of signal events, L int is the integrated luminosity, ǫ is the selection efficiency, and B = B(χ c1,c2 → γJ/ψ) × B(J/ψ → ℓ + ℓ − ) is the branching fraction of intermediate states in the sequential decay, (1 + δ) is the ISR correction factor [29]. The ISR correction factor is calculated with the kkmc program, with the measured √ s-dependent cross section of the reactions e + e − → γχ c1,c2 as input. This procedure is iterated several times until (1+δ)ǫ converges, i.e. the relative difference between the last two iterations is less than 1%. The final measured cross sections σ( √ s) for e + e − → γχ c1 are shown in Fig. 2 and are summarized in Table A in the Appendix. Note that some of the cross sections are negative, thereby seem unphysical. This is caused by the fact that the number of events in the signal region is less than the estimated number of background events from sideband regions, and it can be explained by statistical fluctuations. To study the possible resonances in the e + e − → γχ c1 process, a maximum likelihood fit is performed to the √ s-dependent cross sections. To describe the data, we use two coherent Breit-Wigner (BW) resonances, i.e. the ψ(4040) and ψ(4160), together with a continuum term and the incoherent ψ(3686), ψ(3770)   tail contributions. Since the contribution from ψ(3686) and ψ(3770) is small at √ s > 4 GeV and also lack of data between 3.77 and 4.0 GeV, we do not consider the interference effect from ψ(3686) and ψ(3770). The possible interference effect between continuum and other components is also investigated and we find its contribution is small (and taken as systematic effects). The fit function is thus written as where φ i is the relative phase of the amplitude, and A cont is the continuum amplitude which is parametrized as where f cont and n are the free parameters. BW function is described as where M, Γ tot R and Γ ee R are the mass, full width and electric width of the resonance R, respectively. B R is the branching fraction of R → γχ c1 , and Φ( √ s) is the phase space factor. For the E1 transition of the process ψ(3686) → γχ c1 , we consider an additional factor the E 3 γ [30] and a damping factor [31], according to where the phase space factor is given by Φ( √ s) = Eγ E 0 γ 3/2 [30], and the damping factor as D( [31]. The parameters E γ and E 0 γ are the energy of the E1 photon for the ψ(3686) → γχ c1 decay at c.m. energy √ s and at the ψ(3686) mass, respectively.
In the cross section fit, the likelihood function is defined as where N sig is the number of signal events measured from dataset i, s i is the expected number of signal events for the corresponding dataset, and G represents a Gaussian distribution which describes datasets with high statistics at √ s = 3.773 and 4.178 GeV. N obs j is the number of events observed in the χ c1 mass interval from dataset j, (s + b) j is the expected sum of signal and background events in the same interval, and P represents a Poisson distribution which describes low statistics datasets at other c.m. energies. In the fit PDF, the masses and widths of ψ(3686), ψ(3773), ψ(4040), ψ(4160), and Γ ee · B i [ψ(3686) → γχ c1 ] are fixed to PDG values [12]. The fit result is shown in Fig. 2, and also summarized in Table 1. The significance of ψ(4040), ψ(4160) and the continuum term are estimated to be 3.7σ, 3.3σ and 6.7σ, respectively. Considering the constructive and destructive interferences between ψ(4040) and ψ(4160), there are two solutions with equal good quality from the fit, which has been proved mathematically [32]. A χ 2 -test is used to estimate the fit quality. Due to the low statistics of data at some c.m. energies, we merge the datasets into 17 groups, and the χ 2 -test gives χ 2 /ndf = 10.6/11 = 0.96, where ndf is the number of degree of freedom.   Table 1. Results of the fit to the e + e − → γχc1 cross sections. The unit of e + e − partial width is eV/c 2 and the unit of fcont is eV n /pb. The errors are statistical only.

Parameter
Solution For the e + e − → γχ c2 process, the measured cross sections are shown in the Fig. 3 and summarized in Table B in the Appendix. In the fit PDF, the resonance parameters of ψ(3686), ψ(3770), ψ(4040), and ψ(4160) are also fixed to PDG values. To describe the √ s-dependent cross section, one more resonance is added in the fit function to describe the structure around √ s = 4.39 GeV. According to the fit, the contribution from continuum is not significant in this process (< 1σ). Thus, we construct the fit function as The fit results are shown in Fig. 3 and summarized in Table 2. Same as before, there are four solutions with equal fit quality, due to the interferences between ψ(4040), ψ(4160), and the new resonance [32]. The significance of ψ(4040), ψ(4160) and the resonance near 4.39 GeV are estimated to be 2.0σ, 4.6σ and 5.8σ, respectively. Similarly, we merge the data into 17 groups when performing a χ 2 -test. The χ 2 -test to the fit quality gives χ 2 /ndf = 7.8/9 = 0.87.

A. Event selection
For the e + e − → γχ c0 study, the χ c0 resonance is reconstructed with 2(π + π − ), π + π − K + K − , and K + K − decay modes. Considering the relatively small branching fractions from the χ c0 decay and also the high background levels, only the data samples with L int > 400 pb −1 at √ s > 4.0 GeV are used in this study. The selection criteria of charged tracks and photons are the same as for the e + e − → γχ c1,c2 analysis. The particle identification (PID) of kaons and pions is based on the dE/dx and TOF information, and the particle type with the highest probability is assigned to each track. For photons, the most energetic photon is regarded as the candidate for signal events. A 4C kinematic fit is performed to these three decay modes and χ 2 4c < 25 is required for both χ c0 → 2(π + π − )/K + K − π + π − modes and χ 2 4c < 30 for the χ c0 → K + K − mode.
For the χ c0 → K + K − decay mode, there are Bhabha background events. We require the deposited energy in the EMC over the momentum of a charged track E EMC /p < 0.8 to reject them.
Same as before, |M (γγ extra ) − m(π 0 )| > 10 MeV/c 2 is required to suppress background with π 0 → γγ. Figure 4 shows the M (K + K − π + π − ), M (π + π − π + π − ), and M (K + K − ) invariant mass distributions for the full datasets after imposing the above selection criteria. To obtain the number of χ c0 signal events, an unbinned maximum likelihood fit is performed to the M (K + K − π + π − ), M (π + π − π + π − ) and M (K + K − ) distributions simultaneously at each c.m. energy. The signal yields for three decay modes are constrained according to corresponding reconstruction efficiencies and branching fractions. In the Table 2. Results of the fit to the e + e − → γχc2 cross sections. The unit of the e + e − partial width is eV/c 2 . The errors are statistical only.

Parameter
Solution I Solution II Solution III Solution IV Γ ee B(ψ(3770) → γχc2) fit, the signal PDFs are described with the shapes from simulated signal MC events. The background shapes are described with two 2nd-order polynomial functions for the χ c0 → K + K − π + π − , 2(π + π − ) decay modes, and a 1st-order polynomial function for the K + K − mode. The significance of χ c0 signal is estimated to be less than 2σ at each c.m. energy point. To estimate an upper limit (UL) of the production cross sections, we scan the likelihood curve in the fit and set the 90% C.L. The corresponding UL of the cross section is calculated as where N up is the UL of the number of signal events at 90% C.L., which is obtained by integrating the likelihood curve of the fit (the systematic uncertainty is considered by convolving the likelihood curve with a Gaussian and its standard deviation is set to the systematic uncertainty). (1 + δ v ) is the vacuum polarization factor taken from calculation [33]. The UL of the cross sections at all c.m. energies are shown in Fig. 5 and summarized in Table C in the Appendix .

V. SYSTEMATIC UNCERTAINTY
The systematic uncertainty of the cross section measurements of e + e − → γχ cJ (J = 0, 1, 2) mainly comes from the luminosity measurement, detection efficiency, decay branching fractions, signal extraction, and radiative correction. The luminosity is measured using Bhabha events and the uncertainty is estimated to be 1% [18]. For high momentum leptons, the uncertainty of the tracking efficiency is 1% per track [15]. The uncertainty in the photon reconstruction is 1% per photon, estimated by studying the J/ψ → ρ 0 π 0 decay [34]. The PID efficiency uncertainty for each charged track is taken as 1% [35]. For the systematic uncertainty from the kinematic fit, we correct the track helix parameters in the MC simulation according to the method described in Ref. [36], and the efficiency difference before and after correction is considered as the systematic uncertainty. The uncertainties for the branching fractions of χ c1,c2 → γJ/ψ, and χ c0 → π + π − π + π − , K + K − π + π − , K + K − from the PDG [12] are taken as systematic uncertainties for the cross section measurement. For the e + e − → γχ c0 process, the systematic uncertainties for tracking, PID, photon detection and kinematic fit are the same, and the total systematic error is obtained by weighting each individual one according to the branching fractions and efficiencies of three χ c0 decay modes, by considering the possible correlations between them. For the systematic uncertainty from the background veto requirements, we select the ψ(3686) → γχ c0,c1,c2 control samples, and the selection requirements are exactly the same as the requirements described above. We take the efficiency difference between MC simulation and corresponding control samples as the systematic uncertainties. For the systematic uncertainty from the J/ψ mass window, an e + e − → ηJ/ψ control sample is studied and we take the efficiency difference between MC simulation and control samples as the uncertainties. Table 3. Summary of systematic uncertainties sources on the cross section measurement in %, the "−" indicates that the uncertainty is not applicable. To obtain the systematic uncertainty from signal extraction, we refit the M (γ H J/ψ) by replacing the background shape with 2nd-order polynomial function, varying the fit range (+0.05 GeV/c 2 ), changing the signal shape from signal MC shape to signal MC shape convoluted with a float Gaussian function. The difference with nominal fit results is taken as the systematic uncertainty. For the systematic uncertainties from the ISR correction factor in e + e − → γχ c1,c2 , two sources are considered. First, the difference of (1+δ) * ǫ between the last two iterations is taken as systematic uncertainty, which is 1%. The other sources are the uncertainty of the fit parameters, fit components and damping factor. We sample the fit parameters with Gaussian functions (take the fit results as mean values and the errors as standard deviations) 200 times, then calculate the corresponding (1 + δ) * ǫ values. The standard deviation for (1 + δ) * ǫ of the 200 samplings is taken as a systematic error, which is 1.0% and 1.8% for e + e − → γχ c1 and γχ c2 , respectively. Due to the low significance of ψ(4040)/ψ(4160) in the e + e − → γχ c1 process and ψ(4040) in the e + e − → γχ c2 process, the input cross section line shapes in MC generation are tested by excluding these charmonium states. The differences for (1 + δ) * ǫ are 3.5% and 3.0% for e + e − → γχ c1 and γχ c2 , respectively. For the systematic uncertainty from the damping factor, we change the damping factor from [37](the value of β is also quoted from Ref. [37]), the (1 + δ) * ǫ difference with two different damping factor are 0.9% and 1.3% for e + e − → γχ c1 and γχ c2 , respectively. For the possible interference between continuum and other components in the e + e − → γχ c1 process, the (1 + δ) * ǫ difference with or without considering interference is taken as systematic error, which is 2.3%. For the e + e − → γχ c0 , the difference between a flat line shape or a ψ(3770) line shape is taken as uncertainty. For χ c0 → π + π − π + π − , K + K − π + π − decay modes, the signal MC samples are generated by including all subprocesses. The difference with a pure phase space model is taken as the uncertainty due to the decay model. Table 3 summarizes all the systematic uncertainty sources and their contributions. The total systematic uncertainty is obtained by adding all sources in quadrature.
For the resonance parameters of the structure around 4.39 GeV, the uncertainty of c.m. energies (±0.8 MeV) are common for all data samples, and this uncertainty will propagate directly to the mass measurement. For the uncertainty from the damping factor, we change the damping factor from [31] to e − E 2 γ 8β 2 [37], and the differences are 0.9 MeV and 0.5 MeV for the mass and width, respectively. To estimate the uncertainties from the parameters of ψ(3686), ψ(3770), ψ(4160) and ψ(4040), we randomly sample the four resonances parameters with Gaussian functions (PDG means and errors), and these values are used as input to refit the cross section. The standard deviations of these 1500 fit results are quoted as systematic errors, which are 1.0 MeV/c 2 and 1.6 MeV for the mass and width, respectively. In the cross section fit, the ψ(3770) contribution is added incoherently in the PDF. The possible systematic from the interference effect of ψ(3770) is estimated by considering interference effect in the cross section fit. The differences are 0.8 MeV/c 2 and 0.8 MeV for the mass and width. Assuming all the systematic errors are independent, the total systematic errors are 1.8 MeV/c 2 and 1.9 MeV for mass and width, respectively, by adding all sources in quadrature.

VI. SUMMARY
In summary, using 19.3 fb −1 data at c.m. energies between 3.773 and 4.600 GeV, we observe the e + e − → γχ c1,c2 processes for the first time at √ s = 4.178 GeV.