Study of $e^+e^-\to\gamma\phi J/\psi$ from $\sqrt{s}=4.600$ to $4.951$ GeV

Using data samples with an integrated luminosity of $6.4$~fb$^{-1}$ collected by the BESIII detector operating at the BEPCII storage ring, the process of $e^+e^-\to\gamma\phi J/\psi$ is studied. The processes of $e^+e^-\to\phi\chi_{c1,c2}$, $\chi_{c1,c2}\to\gamma J/\psi$ are observed with a significance of more than $10\sigma$. The $\sqrt{s}$-dependent cross section of $e^+e^- \to \phi\chi_{c1,c2}$ is measured between 4.600 and 4.951~GeV, and evidence of a resonance structure is found for the first time in the $\phi\chi_{c2}$ process. We also search for the processes of $e^+e^-\to\gamma X(4140)$, $\gamma X(4274)$ and $\gamma X(4500)$ via the $\gamma\phi J/\psi$ final state, but no obvious structures are found. The upper limits on the production cross section times the branching fraction for these processes at the 90% confidence level are reported.

At present, the inner structure of these Y -states remains unclear. Experimentally, the e + e − annihilation process is one of the most effective ways to probe the nature of Ystates. The Y (4660) resonance was first observed by the Belle Collaboration in e + e − → π + π − ψ(3686) process via initial-state-radiation (ISR) [5], and subsequently confirmed by the BaBar [6] and BESIII Collaborations [8] in the same process. In the Y (4660) → π + π − ψ(3686) decay, the π + π − system is found to be dominated by a f 0 (980) which has a significant ss component. Recently, the Belle experiment reported the first Y (4626) resonance coupling to the D + s D − s1 (2536) + c.c. meson pair with a significance of 5.9σ [13]. Belle also reported evidence (3.4σ) for a resonance with mass and width consistent with Y (4626) in the D + s D * − s2 (2573) + c.c. process [14]. It is not clear whether Y (4660) and Y (4626) correspond to the same resonance or not. The observation of the Y (4660) state coupling to f 0 (980)ψ(3686) and the Y (4626) state coupling to the charmed-antistrange and anticharmed-strange meson pair may indicate that Y (4660)/Y (4626) have ccss components [15][16][17]. In such a case, the Y (4660)/Y (4626) may also decay to the final states of φχ c1 or φχ c2 . The BESIII experiment has measured the cross section of e + e − → φχ c1,c2 at 4.600 GeV [18], and significant χ c1,c2 signals were found. With the data taken at center-ofmass (c.m.) energies up to √ s = 4.951 GeV at BESIII, which fully covers the Y (4626) and Y (4660) mass region, we are able to measure the cross section line shape of e + e − → φχ c1,c2 . The measurements may shed light on the inner structure of the Y (4660)/Y (4626) states and help us understand their nature.
In addition to the Y -states, the non-vector X-states in the φJ/ψ system also attract much interest. A narrow (Γ = 11.7 MeV) near-threshold peak around 4143 MeV/c 2 in the φJ/ψ mass spectrum was first reported by the CDF Collaboration in the B + → φJ/ψK + process with 3.8σ evidence (labeled as X(4140)) [19]. From the potential model, charmonium states within this mass region are expected to have much larger widths due to the open charm decay channels [11]. The X(4140) is therefore suggested to be a candidate of an exotic state. An updated analysis by the CDF Collaboration in 2011 [20] not only confirmed the existence of X(4140) with a 5σ observation, but also reported evidence (3.1σ) of a new narrow peak near 4274 MeV/c 2 in the φJ/ψ spectrum. Subsequent measurements were also carried out by the Belle [21], LHCb [22][23][24][25], CMS [26], D0 [27,28], BaBar [29] and BESIII experiments [18,30]. Belle [21], BaBar [29] and BESIII [18,30] found no evidence for the X(4140).
In this article, we present a study of the e + e − → γφJ/ψ process with 6.4 fb −1 of data taken at e + e − c.m. energies from 4.600 to 4.951 GeV [41][42][43]. The √ s-dependent cross section of e + e − → φχ c1,c2 is measured and possible vector resonances are investigated. The χ c1,c2 resonances are reconstructed with the γJ/ψ and J/ψ → ℓ + ℓ − (ℓ = e, µ) decays. We also search for the possible X-state in the e + e − → γX, X → φJ/ψ process. To increase the number of candidates, both φ → K + K − and φ → K 0 S K 0 L modes are used to reconstruct φ.

BESIII detector and MC sample
The BESIII detector [44] records symmetric e + e − collisions provided by the BEPCII storage ring [45], which operates with a peak luminosity of 1×10 33 cm −2 s −1 at e + e − center-of-mass energy 3.77 GeV. BESIII has collected large data samples between 2.0 and 4.951 GeV [46]. 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 chamber (MUC) system 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. The time resolution in the TOF barrel region is 68 ps, while that in the end cap region was 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 [47]. All the data sets with √ s > 4.600 GeV are taken with the new end cap TOF system.
Simulated samples produced with a geant4-based [48] Monte Carlo (MC) software, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds. The signal MC samples of e + e − → φχ c1,c2 → γφJ/ψ and e + e − → γX → γφJ/ψ are simulated at each c.m. energy point corresponding to the luminosity of data, with φ → K + K − , K 0 S K 0 L and J/ψ → e + e − , µ + µ − being simulated according to the branching fractions taken from the Particle Data Group (PDG) [10]. The inclusive MC sample includes the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes simulated with kkmc [49]. The simulation models the beam energy spread and ISR in the e + e − annihilations with the generator kkmc [49]. The known decay modes of charmed hadrons are modelled with evtgen [50] using branching fractions taken from the PDG [10], and the remaining unknown decays are modelled with lundcharm [51]. Final state radiation (FSR) from charged final state particles is incorporated using photos [52].

Event Selection
For candidate events of interest, the φ meson is reconstructed from S is reconstructed from π + π − and the K 0 L is missing due to the low detection efficiency with the BESIII detector. The χ c1,c2 is reconstructed from γJ/ψ, and the J/ψ is reconstructed from the lepton pairs e + e − or µ + µ − . The following event selection criteria are applied to both data and MC samples.
Charged tracks detected in the MDC are required to be within a polar angle (θ) range of |cosθ| < 0.93 (the coverage of the MDC), where θ is defined with respect to the z-axis, which is the symmetry axis of the MDC. For charged tracks not used for K 0 S reconstruction, the distance of closest approach to the interaction point (IP) must be less than 10 cm along the z-axis, |V z | < 10 cm, and less than 1 cm in the transverse plane, |V xy | < 1 cm, while those for K 0 S reconstruction, only a loose requirement of |V z | < 20 cm is applied. Photon candidates are identified using showers in the EMC. The deposited energy of each shower must be greater than 25 MeV in the barrel region (| cos θ| < 0.80) and greater than 50 MeV in the end cap region (0.86 < | cos θ| < 0.92). To exclude the showers that originate from charged tracks, the angle between the position of each shower in the EMC and the closest extrapolated charged track must be greater than 10 degrees. To suppress the electronic noise and the showers unrelated to the event, the difference between the EMC time and the event start time is required to be within [0, 700] ns.
For each event, the lepton pair (ℓ + ℓ − ) from J/ψ decays and the kaons from φ decays can be effectively distinguished by their momenta in the lab-frame. The tracks with momentum larger than 1 GeV/c are assigned as leptons, and the amount of deposited energy in the EMC is further used to separate the muons from electrons. For both muon candidates, the deposited energy in the EMC is required to be less than 0.4 GeV, while it is required to be greater than 1.0 GeV for electrons. For the tracks with momentum less than 1 GeV/c, particle identification (PID), which combines measurements of the energy deposited in the MDC (dE/dx) and the flight time in the TOF to form likelihoods L(h) (h = K, π) for each hadron h hypothesis, is used. Tracks are identified as kaons when the kaon hypothesis has a larger likelihood than the pion hypothesis (L(K) > L(π) and L(K) > 0).

3-track events with
For the φ → K + K − channel, one of the kaons could be missing due to an inefficiency. Together with the lepton pair from the J/ψ decay, there are three charged particles remaining in each signal event (referred to as the 3-track events). Two of the charged tracks are assigned as the lepton pair and the third as a kaon. The PID likelihood of the kaon is required to satisfy L(K) > L(π) and L(K) > 0, and at least one photon candidate is also required.
To improve the resolution and suppress background, a one-constraint (1C) kinematic fit is applied to the 3-track event by constraining the mass of the missing particle to the kaon nominal mass (M (K ∓ miss ) = (P e + e − − P γK ± ℓ + ℓ − ) 2 ) inferred from the four momentum conservation. For the events with multiple photons in the final state, the combination of γK ± K ∓ miss ℓ + ℓ − with the smallest χ 2 from the kinematic fit is retained, and χ 2 < 20 is required.
To reduce the π misidentification background in the µ + µ − channel, the MUC is used to identify muons. At least one of muon candidate should have a hit depth > 30 cm in the MUC. To veto the radiative Bhabha background in J/ψ → e + e − events, the polar angle of e + is required to satisfy cos(θ e + ) < 0.85.

4-track events with
For a candidate event with K + K − ℓ + ℓ − detected (referred to as 4-track events), the photon candidate is always ignored and not required to be detected in order to improve the efficiency. At least four charged tracks are required, two of which are assigned as the lepton pair and the remaining tracks as kaons. Both kaons are required to be identified. A similar 1C kinematic fit is performed by constraining the mass of the missing particle to be a photon inferred from the four momentum conservation, i.e. M (γ miss ) = (P e + e − − P K The kinematic fit χ 2 is required to be χ 2 < 35. The same MUC requirement as for 3-track events is applied to the muon candidates to suppress pion background.

Events with
candidate has a long lifetime and is not detected. We require at least four charged tracks to be detected in each event, two of which are assigned as the lepton pair and the remaining charged tracks are pions. The K 0 S candidates are reconstructed from two oppositely charged pions satisfying |V z | < 20 cm. There is no PID requirement for the charged pions, and they are constrained to originate from a common secondary decay vertex. The decay length of the K 0 S candidate is required to be greater than twice the vertex resolution away from the IP to suppress the non-K 0 S background. After the vertex fit, we set a mass window of 0.490 GeV/c 2 < M (π + π − ) < 0.505 GeV/c 2 for the K 0 S candidate (the mass resolution is 5 MeV/c 2 ). At least one good photon candidate is also required in each event.
A 1C kinematic fit is applied to each event, with the mass of the missing particle constrained to the K 0 L nominal mass inferred from the four momentum conservation, i.e.

Cross section measurement
Based on the event selection, the χ c1 and χ c2 signals are observed from both the φ → K + K − and K 0 S K 0 L modes. To determine the signal yields, an unbinned maximum likelihood fit is performed to the M (γJ/ψ) distribution in the φ → K + K − and K 0 S K 0 L modes simultaneously. In the fit at each c.m. energy, the signal probability-density-function (PDF) is described by a MC-simulated shape convolved with a Gaussian function, which models the resolution difference between data and MC simulation. The MC-simulated shape is a weighted sum of the simulations at each c.m. energy, which has already taken into account the c.m. energy and decay modes dependence for the resolution. The Gaussian parameters are determined from the fit to the full dataset which has higher statistics. A linear function is used to describe the background. The two modes share the same φχ c1,c2 production cross section at the same c.m. energy. The selection efficiencies and branching fractions of the φ → K + K − /K 0 S K 0 L modes at each c.m. energy are included in the fit. Figure 4 shows the fit result for the full dataset from √ s = 4.600 GeV to 4.951 GeV, and the corresponding plots at each individual c.m. energy are shown in Figure 9 of Appendix A. The statistical significance is estimated by comparing the fit likelihoods with and without the χ c1,c2 signal. In addition to the nominal fit, the fits by changing the background shape and the fit range have also been performed. In all the cases, the significance of the χ c1,c2 is found to be greater than 10σ, by comparing the difference of log-likelihoods ∆(−2 ln L) = 137(131) for the χ c1 (χ c2 ) and taking into account the change of the number of degrees of freedom (∆d.o.f. = 5).  Figure 4. Sum of the simultaneous fits to M (γJ/ψ) distribution for the full data sets. Dots with error bars are the data, the solid curves are the fit results, the red dashed, blue dotted, and green dash-dotted lines are the χ c1 , χ c2 and the background shape, respectively.
The Born cross section of e + e − → φχ c1,c2 at c.m. energy √ s is calculated with where N fit is the number of fitted events for the φχ c1,c2 , which is equal to the number of the φχ c1,c2 events in data divided by the efficiency and branching fraction of φ, L int is the integrated luminosity, (1 + δ) is the ISR correction factor obtained from kkmc, 1 |1−Π| 2 is the vacuum polarization factor [53], and B is the product of the branching fraction for χ c1,c2 → γJ/ψ and J/ψ → ℓ + ℓ − . The Born cross sections of e + e − → φχ c1 and e + e − → φχ c2 at each c.m. energy are listed in Tables 1 and 2, respectively. In case the signal significance is less than 3σ, an upper limit of the Born cross section (σ U.L. ) at the 90% confidence level (C.L.) is also reported. The upper limits of φχ c1 and φχ c2 yields are estimated via a Bayesian approach [10]. A likelihood scan L(n) is performed with various assumptions for the number of signal events (n) in the fit. The systematic uncertainty is also considered by smearing the likelihood distribution with a Gaussian function with width equal to the systematic uncertainty. The upper limit of N U.L. at the 90% C.L. corresponds to The detection efficiencies of φχ c1,c2 events depend on the e + e − c.m. energy. With the increasing c.m. energy, charged kaons have higher momentum and are thus much more efficient to be detected, while for the K 0 S K 0 L channel, due to more ISR events the reconstruction efficiency whereas decreases (the π + π − from K 0 S decay already have sufficient momentum to be detected and are not sensitive to c.m. energy).
L for the 3-track events of the φ → K + K − mode, 4-track events of the φ → K + K − mode and the events of the φ → K 0 S K 0 L mode, respectively, the product of radiative correction factor and vacuum polarization factor 1+δ |1−Π| 2 and the number of fitted events N fit (also the corresponding upper limit N U.L. at the 90% C.L. in case the signal significance is less than 3σ). The first uncertainty is statistical and the second is systematic.  at the 90% C.L. in case the signal significance is less than 3σ. The table also includes integrated luminosity L int , detection efficiency ǫ 3 K + K − , ǫ 4 K + K − and ǫ K 0 S K 0 L for the 3-track events of the φ → K + K − mode, 4-track events of the φ → K + K − mode and the events of the φ → K 0 S K 0 L mode, respectively, the product of radiative correction factor and vacuum polarization factor 1+δ |1−Π| 2 and the number of fitted events N fit (also the corresponding upper limit N U.L. at the 90% C.L. in case the signal significance is less than 3σ). The first uncertainty is statistical and the second is systematic.
To investigate the √ s-dependent cross section line shape of e + e − → φχ c1,c2 , a maximum likelihood fit is performed to the dressed cross section (σ B ( √ s) 1 |1−Π| 2 ). Due to the small numbers of events at each single c.m. energy, the likelihood function is constructed as where P represents a Poisson distribution, N obs i , N exp i and N bkg i are the number of observed events, the number of expected χ c1,c2 signal events and the background events in the χ c1,c2 signal region for the i-th dataset, respectively. Here in the fit, only statistical uncertainties are considered.
For the e + e − → φχ c1 process, a continuum amplitude is used to fit the cross section, where f cont and n are free parameters in the fit. We also use a phase space (PHSP) shape corrected continuum amplitude A cont ( √ s) Φ( √ s) to fit the cross section, where Φ( √ s) is the two-body PHSP factor. Figure 5 shows the fit results with both models, and the numerical results are listed in Table 3. We also fit the cross section data with a Breit-Wigner (BW) function and the coherent sum of a BW and a continuum amplitude, and no significant resonance structures are found.   For the e + e − → φχ c2 process, there is a possible resonance structure around 4.7 GeV in the cross section line shape as shown in Fig. 6, which is fitted with a BW function: where M , Γ tot and Γ e + e − are the mass, full width, and electric width of the potential resonance Y , respectively, and B(Y → φχ c2 ) is the branching fraction of Y → φχ c2 . for the resonance. A χ 2 test method is used to estimate the fit quality, which gives χ 2 /d.o.f. = 15.9/10. The significance for the resonance hypothesis over the continuum hypothesis is estimated to be 3.1σ, by comparing the difference of log-likelihoods ∆(−2 ln L) = 10.0 and taking into account the change of number of degree of freedom (∆d.o.f. = 1). Here the continuum hypothesis follows A cont ( √ s) Φ( √ s). The fit result for the continuum hypothesis is shown in Figure 6 (a) (dash-dotted line) and listed in Table 4 (last column). The potential resonance (solid line in Figure 6 (a)) is found to be consistent with the Y (4660) reported in e + e − → π + π − ψ(3686) [5][6][7][8]. Next, we fit the e + e − → φχ c2 cross section with the fixed mass and width of the Y (4660) [8]. Two fit models are considered: one is the single BW model, which gives Γ e + e − B[Y (4660) → φχ c2 ] = 1.0 ± 0.1 eV with a fit quality χ 2 /d.o.f. = 21.5/12 (the dashed line in Figure 6 Figure 6 (a)).
Since the fit quality with the fixed Y (4660) is close to the one with a single free BW model (χ 2 /d.o.f. = 15.9/10), we cannot distinguish between these two models.
To improve the fit quality, the fit model is parameterized as the coherent sum of a BW resonance and a possible continuum term (BW + A cont e iφ ). The fit result is shown in Figure 6 Table 4. The significance for the coherent sum of a BW and continuum model (BW + A cont e iφ ) over the single BW model is estimated to be 2.3σ. Thus, we are not able to distinguish these two models based on the current data.
Parameter Since no obvious structures are observed in the φχ c1 mode, the upper limit of Γ e + e − B(Y → φχ c1 ) is also determined for the possible structures observed in the φχ c2 mode. A similar method by scanning the Γ e + e − B(Y → φχ c1 ) dependent likelihood distribution is used, and the results at 90% C.L. are listed in Table 5.

Systematic uncertainty for cross section measurement
The sources of systematic uncertainties in the cross section measurement of e + e − → φχ c1,c2 include the luminosity measurement, tracking efficiency, PID efficiency, K 0 S reconstruction, photon reconstruction, kinematic fit, radiative correction, MC model, MUC response, branching ratios, and the fit.
The uncertainty of the integrated luminosity measurement is 0.6% by analyzing the large angle Bhabha events at BESIII [42]. The uncertainty of the tracking efficiency for high momentum leptons is 1% per track, and thus 2% by adding both leptons linearly [54] since we require both leptons detected. For the φ → K + K − mode, both one kaon events and two kaon events are reconstructed. Assuming p (q) is the corresponding tracking efficiency for a single kaon from data (MC), the efficiency to reconstruct both one and two kaon candidates is 2p(1 − p) for data (MC). Considering p ≈ 85% and the tracking efficiency uncertainty p/q − 1 = 1% at BESIII, the uncertainty due to the detection of both one and two kaon candidates for the tracking efficiency can be calculated as 1 − 1−(1−p) 2 1−(1−q) 2 , which is negligible. The same calculation can be applied to the kaon PID uncertainty, which is also negligible. For the φ → K 0 S K 0 L mode, the uncertainty of tracking efficiency is 1% per pion. The uncertainty of K 0 S reconstruction is estimated to be 1.2% by studying the J/ψ → K * (892) ± K ∓ → K 0 S π ± K ∓ and J/ψ → φK 0 S K ∓ π ± control samples [55]. The uncertainty from photon reconstruction is estimated to be 1% per photon by studying the J/ψ → ρ 0 π 0 events [56].
The systematic uncertainty associated with kinematic fitting is estimated by comparing the efficiency difference with or without the helix parameters correction in MC simulations [57]. The radiative correction factor and efficiency depend on the input cross section line shape in kkmc. Using different cross section line-shapes as studied in section 3.2, the difference in (1 + δ)ǫ between different models is taken as the systematic uncertainty. In the signal MC simulation, a phase space model is used. To estimate the uncertainty due to the MC model, the angular distribution of e + e − → φχ c1,c2 is modelled by a 1 ± cos 2 θ distribution, and the efficiency difference is taken as the systematic uncertainty.
The uncertainty from the MUC response is studied with a control sample of e + e − → µ + µ − events. The difference in efficiency between the data and MC simulation due to the requirement of µ hit depth in the MUC is taken as the systematic uncertainty. The uncertainties of branching fractions of the intermediate states are taken from the PDG [10]. The uncertainties related to the fit are investigated by changing the fit range and changing the background shape from a free 1st-order polynomial to a fixed flat shape with the number of events estimated from φ and J/ψ sidebands. The largest difference in signal yields is taken as the systematic uncertainty.
In section 3, three data samples, which are the 3-track events, 4-track events with φ → K + K − and the events with φ → K 0 S K 0 L , are reconstructed. A source of systematic uncertainty can contribute differently to the three data samples. To propagate the systematic uncertainty to the cross section, we take the weighted average of the systematic uncertainties in the three data samples, which follows where σ tot is the average systematic uncertainty to the cross section as listed in Table 6, ω i and σ i are the weight and systematic uncertainty for ith data sample, ǫ i and B i are the efficiency and branching ratio of φ for the ith data sample, ρ ij is the correlation parameter between the ith and jth data samples, and ρ ij = 1 if the systematic uncertainty is correlated between the ith and jth data samples, otherwise ρ ij = 0. Assuming all these sources are independent, the total systematic uncertainty in the cross section measurement is obtained by adding them in quadrature. Table 6 summarizes all the systematic sources and their contributions at 4.68 GeV, and the systematic uncertainties at other energy points are listed in Tables 13 and 14

Systematic uncertainties for the resonance parameters
The systematic uncertainties for the resonance parameters mainly come from the absolute c.m. energy calibration, the parameterization of the BW function, and the cross section measurement. The c.m. energies of the data sets used in this work are measured with Λ c events, with an uncertainty of ±0.6 MeV [42,43]. This common uncertainty for all the data samples could shift the cross section line-shape globally, and is thus the systematic uncertainty to the mass of the resonance.
In the fit to the cross section of e + e − → φχ c2 (Figure 6), a constant full width BW function is employed. We also use an alternative BW function, where the constant width is replaced by an energy dependent width Γ( Here Γ 0 is the full width at √ s = M . The difference in the resonance parameters between the two BWs is taken as the systematic uncertainty. In the fit to the cross section of e + e − → φχ c2 (Figure 6 (b)), a continuum amplitude (Equation (3.3)) is used to describe the non-resonance contribution. We also use a PHSP corrected continuum amplitude (A cont √ Φ) in the fit. The difference in the resonance parameters is taken as the systematic uncertainty.
The uncertainty from the cross section measurement can be divided into two parts, one is the correlated systematic uncertainty for all the energy points, including tracking, photon reconstruction, K 0 S reconstruction, luminosity, branching fraction, muon hit depth, background shape, and fit range. They are propagated to Γ e + e − B(Y → φχ c2 ) directly. The other is the uncorrelated systematic uncertainty, which is dominated by the radiation correction according to the previous section. This uncertainty can be considered in the fit to the cross section. The two types of uncertainties are added in quadrature assuming they are independent. Tables 7 and 8 summarize the sources of systematic uncertainty for the resonance parameters and their contributions, and the total systematic uncertainty is obtained by adding them in quadrature.

Source
Mass 4 Study of e + e − → γX with X → φJ/ψ The process of e + e − → γX → γφJ/ψ shares the same final states as that of e + e − → φχ c1,c2 , thus the same event selection criteria are applied to the e + e − → γX process. The M (φJ/ψ) invariant mass distribution, shown in Figure 7, is well described by the φχ c1,c2 events, together with the non-γφJ/ψ background events estimated from the φ-J/ψ 2-dimensional sidebands (B1/2 + B2/2 − B3/4 as exhibited in Figures 1 to 3). No other structure is observed in the M (φJ/ψ) mass distribution.
-18 -  Figure 7. The invariant mass distribution of M (φJ/ψ) in the φ → K + K − and φ → K 0 S K 0 L modes. Dots with error bars are the full data, the red dashed and blue dotted histograms are from φχ c1 and φχ c2 MC, which have been normalized to the data, the black solid histograms are the sum of φχ c1 and φχ c2 , and the green filled histograms are the φ − J/ψ 2-dimensional sideband.

Upper limit of e + e − → γX cross section
The product of Born cross section of e + e − → γX and the branching fraction of X → φJ/ψ is calculated by where N fit γX is the number of fitted events for γX, which is equal to the number of γX events in data divided by the efficiency and branching fraction of φ, L int is the integrated luminosity, 1 + δ is the ISR correction factor, 1 |1−Π| 2 is the vacuum polarization factor, and B is the branching fraction of J/ψ → ℓ + ℓ − .
Since no significant structures are observed, we determine the upper limit of the production cross section for e + e − → γX → γφJ/ψ using the same method as described in section 3.2. An unbinned maximum likelihood fit is performed to the M (φJ/ψ) distribution simultaneously for the φ → K + K − and K 0 S K 0 L modes. In the fit, the signal PDF is described by MC-simulated shapes, where the mass and width of X are fixed to LHCb's measurements [25]. The background is composed of φχ c1,c2 and a smooth polynomial shape (including both the non-γφJ/ψ and the continuum γφJ/ψ contribution). The φχ c1,c2 background shapes are from the MC simulation, and their yields are normalized to the cross section measurement described in section 3.2. The contribution for the sum of non-γφJ/ψ and continuum γφJ/ψ backgrounds is free. The selection efficiencies and branching fractions of φ → K + K − /K 0 S K 0 L modes are also included in the fit procedure. Figure 8 shows the upper limit of the Born cross section at the 90% C.L. for e + e − → γX → γφJ/ψ at each c.m. energy, and the numerical results are listed in Tables 9 to 11.   Table 9. The upper limit of Born cross section at 90% C.L. σ U.L. B(X → φJ/ψ) for e + e − → γX(4140) at each c.m. energy √ s. The table also includes integrated luminosity L int , detection efficiency ǫ 3 K + K − , ǫ 4 K + K − and ǫ K 0 S K 0 L for the 3-track events in the φ → K + K − mode, 4-track events in the φ → K + K − mode and the events in the φ → K 0 S K 0 L mode, respectively, the product of radiative correction factor and vacuum polarization factor 1+δ |1−Π| 2 , and the 90% C.L. upper limit of the number of fitted events for γX(4140) N U.L. γX(4140) .
L for the 3-track events in the φ → K + K − mode, 4-track events in the φ → K + K − mode and the events in the φ → K 0 S K 0 L mode, respectively, the product of radiative correction factor and vacuum polarization factor 1+δ |1−Π| 2 and the 90% C.L. upper limit of the number of fitted events for γX(4274) N U.L.
L for the 3-track events in the φ → K + K − mode, 4-track events in the φ → K + K − mode and the events in the φ → K 0 S K 0 L mode, respectively, the product of radiative correction factor and vacuum polarization factor 1+δ |1−Π| 2 and the 90% C.L. upper limit of the number of fitted events for γX(4500) N U.L. γX(4500) .

Systematic uncertainty
Since the same selection criteria have been applied to the e + e − → φχ c1,c2 and e + e − → γX processes, they share most of the systematic uncertainties, such as the tracking efficiency, PID efficiency etc. (cf. section 3.3), and their contributions are listed in Table 12. The systematic uncertainties specifically for the e + e − → γX process are described below.
The uncertainty due to the signal shape is considered by varying the mass and width of X states within ±1σ, and changing the signal shape to a MC shape convolved with a 2 MeV Gaussian resolution function [58]. For the uncertainty due to background, the number of φχ c1,c2 background events is varied within ±1σ, and the smooth polynomial background is studied by varying the order of the polynomial or replacing it with a shape estimated from the sideband data in the fit. The uncertainty associated with the fit range is determined by varying the fit range within ±10 MeV. By taking these sources into consideration in the fit, the most conservative upper limit for e + e − → γX is reported.
We search for potential vector Y -states in the cross section line shape of e + e − → φχ c1,c2 , which might contain ccss components in their internal structure. For the e + e − → φχ c1 process, we find no obvious structure in the cross section line shape, and a continuum amplitude can well describe it. For the e + e − → φχ c2 process, there is an enhancement in the cross section line shape. A fit to the cross section with a single BW resonance gives M = (4672.8±10.8±3.9) MeV/c 2 and Γ = (93.2±19.8±9.4) MeV for the mass and width of the structure. The significance of the resonance hypothesis over non-resonance hypothesis is estimated to be 3.1σ. The mass and width of the resonance are consistent with the Y (4660) reported in e + e − → π + π − ψ(3686) [5][6][7][8]. An alternative fit to the cross section with the coherent sum of a BW and a continuum amplitude gives M = (4701.8 ± 10.9 ± 2.7) MeV/c 2 and Γ = (30.5 ± 22.3 ± 14.6) MeV for the mass and width of the structure, which has a higher mass and narrower width. The significance for the resonance hypothesis in this model is estimated to be 3.6σ. However, within the current uncertainties, we are not able to distinguish whether it is the same structure as the Y (4660), and the significance for the second fit over the first one is only 2.3σ. This is the first evident structure observed in the φχ c2 system. We also search for a possible X-state in the φJ/ψ system through the radiative process e + e − → γX → γφJ/ψ. The φJ/ψ spectrum can be well described by the φχ c1,c2 and background events, and no other structure is evident in the M (φJ/ψ) mass distribution. The X(4140), X(4274) and X(4500) resonances reported by the LHCb Collaboration [25] are not observed, and the upper limits on the Born cross sections for e + e − → γX(4140), γX(4274), γX(4500) → γφJ/ψ at the 90% C.L. are determined.