Cross section measurements of $e^{+}e^{-} \to K^{+}K^{-}K^{+}K^{-} $ and $ \phi K^{+}K^{-}$ at center-of-mass energies from 2.10 to 3.08 GeV

We measure the Born cross sections of the process $e^{+}e^{-} \to K^{+}K^{-}K^{+}K^{-}$ at center-of-mass (c.m.) energies, $\sqrt{s}$, between 2.100 and 3.080 GeV. The data were collected using the BESIII detector at the BEPCII collider. An enhancement at $\sqrt{s}= 2.232$ GeV is observed, very close to the $e^{+}e^{-} \to \Lambda \overline{\Lambda}$ production threshold. A similar enhancement at the same c.m. energy is observed in the $e^{+}e^{-} \to \phi K^{+}K^{-}$ cross section. The energy dependence of the $K^{+}K^{-}K^{+}K^{-}$ and $\phi K^{+}K^{-}$ cross sections differs significantly from that of $e^{+}e^{-} \to \phi \pi^{+}\pi^{-}$.

We measure the Born cross sections of the process e + e − → K + K − K + K − at center-of-mass (c.m.) energies, √ s, between 2.100 and 3.080 GeV. The data were collected using the BESIII detector at the BEPCII collider. An enhancement at √ s = 2.232 GeV is observed, very close to the e + e − → ΛΛ production threshold. A similar enhancement at the same c.m. energy is observed in the e + e − → φK + K − cross section. The energy dependence of the K + K − K + K − and φK + K − cross sections differs significantly from that of e + e − → φπ + π − .

I. INTRODUCTION
The φ(2170) resonance, denoted previously as Y (2175), was first observed by BABAR in the process e + e − → φf 0 (980) → φππ [1] via initial-state radia-tion (ISR) and was confirmed by Belle [2]. BES [3] and BESIII [4,5] also observed the φ(2170) in the φf 0 (980) invariant-mass spectrum. The discovery of ss bound states is of interest for the understanding of the strangeonium spectrum, which is less well understood than for example the hidden-charm states (cc). The CLEO Collaboration found the first evidence for Y (4260) → K + K − J/ψ [6] above the DD-production threshold. A similar process, e + e − → φK + K − , potentially allows the study of strangeonium-like vector states above the KK-production threshold.
Many theoretical interpretations have been proposed for the φ(2170), such as a ssg hybrid [7], a 2 3 D 1 ss state [8], a tetraquark state [9,10], a ΛΛ bound state [11,12], or a three-meson system φK + K − [13]. The 1 −− ssg hybrid can decay to φππ, with a cascade (ssg → (ss)(gg) → φππ) [14], whereby ssg → φf 0 (980) may make a significant contribution. However, none of the theoretical models has so far been able to describe all experimental observations in all aspects. Searching for new decay modes and measuring the line shapes of their production cross sections will be very helpful for interpreting the internal structure of the φ(2170) resonance.
The BABAR Collaboration measured the e + e − → K + K − K + K − cross sections and observed an enhancement around 2.3 GeV [15,16]. In addition, the BES Collaboration observed the f 0 (980), f ′ 2 (1525) and f 0 (1790) in the invariant-mass distribution of K + K − pairs in events in which the other K + K − pair has an invariant mass close to the nominal φ mass [17]. An enhancement at √ s = 2.175 GeV was seen in the line shape of the process e + e − → φf 0 (980) [16], but due to poor statistics, no strong conclusion could be drawn from the data. Torres et al. have performed a Faddeev calculation for the three-meson system φK + K − and obtained a peak around 2.150 GeV/c 2 [13]. These observations stimulate experimentalists to study the energy dependence for the production of the φK + K − and K + K − K + K − final states.
Using a data sample corresponding to an integrated luminosity of 650 pb −1 collected at center-of-mass (c.m.) energies from 2.0 GeV to 3.08 GeV [18], we present in this paper the results of a study of the reaction e + e − → K + K − K + K − and its dominant intermediate process e + e − → φK + K − .

II. DETECTOR AND DATA SAMPLES
The BESIII detector is a magnetic spectrometer [19] located at the Beijing Electron Positron Collider (BEPCII) [20]. The cylindrical core of the BE-SIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorime-ter (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 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 optimization of event-selection criteria, the determination of detection efficiencies and the estimates of potential backgrounds are performed based on Monte Carlo (MC) simulations taking the various aspects of the experimental setup into account. The geant4-based [21] MC simulation software, which includes the geometric and material description of the BESIII detector, the detector response and digitization models, and the detector running conditions and performances, is used to generate the MC samples.
For the background study, the e + e − → qq process is simulated by the MC event generator conexc [22], while the decays are generated by evtgen [23,24] for known decay modes with branching fractions set to Particle Data Group (PDG) world-average values [25] and by luarlw [26] for the remaining unknown decays. MC samples of e + e − → e + e − and µ + µ − processes are generated by babayaga 3.5 [27]. The signal MC samples from the phase-space models (PHSP) of e + e − → K + K − K + K − and e + e − → φK + K − are generated at c.m. energies corresponding to the experimental values, where the line shape of the production cross section of the two processes is taken from the BABAR experiment [16] and the signal detection efficiency is obtained by weighting the MC-generated PHSP sample to data according to the observed invariant-mass distribution.

III. EVENT SELECTION AND BACKGROUND ANALYSIS
A. e + e − → K + K − K + K − Candidate events are required to have three or four charged tracks. Charged tracks are reconstructed from hits in the MDC within the polar angle range |cosθ| < 0.93 and are required to pass the interaction point within 10 cm along the beam direction and within 1 cm in the plane perpendicular to the beam. For each charged track, the TOF and the dE/dx information are combined to form particle identification (PID) confidence levels (C.L.) for the π, K, and p hypotheses. The particle type with the highest C.L. is assigned to each track. At least three kaons are required to be identified. The primary vertex of the event is reconstructed by three kaons. For events with four identified kaons, the combination with the smallest chi-square of the vertex fit is retained. Figure 1 shows the momentum distribution of the three identified kaons for √ s = 2.125 GeV after applying the above-mentioned selection criteria. The peak on the right-side of the spectrum stems from reducible QED background, dominated by the processes e + e − → e + e − and e + e − → µ + µ − . To suppress this background, the momenta of the identified particles are required to be less than 80% of the mean momentum of the colliding beams (p beam ). The black dots with error bars are data, the dashed (red) histogram is from e + e − → qq, the solid (green) histogram is from e + e − → e + e − , the hatched (black) histogram is from e + e − → µ + µ − , and the dotted (blue) histogram is the sum of all MC samples.
The selection criteria for three or four kaons are the same as described in the previous subsection. In addition to the primary-vertex fit of the three kaons, a one-constraint (1C) kinematic fit is performed under the hypothesis that the KK + K − missing mass corresponds to the kaon mass. For events with four reconstructed and identified kaons, the combination with the smallest chi-square of the 1C kinematic fit (χ 2 1C (K + K − KK miss )) is retained and required to be less than 20. In the following, the K miss momentum is that obtained from the 1C kinematic fit and is used in invariant-mass calculations.
The open histogram in Fig. 2 shows the invariantmass distribution for all K + K − pairs for the selected K + K − K + K − events (four entries per event) for data taken at √ s = 3.080 GeV. The hatched histogram in the same figure corresponds to the distribution of the pair with a mass closest to the nominal φ mass. A promi-nent peak near the φ mass is seen in both histograms and indicates that the φ K + K − channel dominates the

IV. SIGNAL YIELDS
The signal yields of e + e − → K + K − K + K − are obtained from unbinned maximum-likelihood fits to the The signal is described by the line shape obtained from the MC simulation convolved with a Gaussian function, where the Gaussian function describes the difference in resolution between data and MC simulation. The background shape is parametrized by a second-order Chebyshev polynomial function. The parameters of the Gaussian function and the Chebyshev polynomial function are left free in the fit. The corresponding fit result for data taken at √ s = 3.080 GeV is shown in Fig. 3.
To determine the signal yields of the e + e − → φK + K − process, an unbinned maximum-likelihood fit is performed to the M (K + K − ) spectra. The probability density function of the M (K + K − ) spectra for the φ is obtained from a P -wave Breit-Wigner function convolved with a Gaussian function that accounts for the detector resolution. The P -wave Breit-Wigner function is defined as where m 0 is the nominal φ mass as specified by the PDG, p is the momentum of the kaon in the rest frame of the K + K − system, p ′ is the momentum of the kaon at the nominal mass of the φ, and Γ 0 is the width of the φ. The angular momentum (ℓ) is assumed to equal one, which is the lowest allowed given the parent and daughter spins, B(p) is the Blatt-Weisskopf form factor, and R is the radius of the centrifugal barrier, whose value is taken to be 3 GeV/c −1 [28].
The background shape is described by an ARGUS function [29]. The parameters of the Gaussian function and the ARGUS function are left free in the fit. The corresponding fit result for data taken at √ s = 3.080 GeV is shown in Fig. 4.
The same event selection criteria and fit procedure are applied to the other 19 data samples taken at different c.m. energies. The number of events for these samples are listed in Tables I and II.
The detection efficiency is obtained by MC simulations of the φK + K − channel using PHSP. It is found that data deviate strongly from the PHSP MC distributions, as demonstrated by the histograms in  which show the non-φ pair K + K − invariant-mass distributions. Here, φ candidates are selected in the signal region and background from the side-band region shown in Fig. 4. The background in Fig. 5 is the distribution of the invariant mass of the remaining pair in the side-band event, and the data points are the invariant mass of the remaining pair of the φ candidates minus the background. To obtain a more accurate detection efficiency, the MC-generated events are weighted according to the observed K + K − (non-φ pair) invariant-mass distribution, where the weight factor is the ratio of the K + K − mass distribution between data and PHSP MC. The weighted PHSP MC distribution is consistent with the background-subtracted data, as shown by the solid histogram in Fig. 5. The detection efficiencies determined by using the weighted MC data and by using the φK + K − PHSP MC data do not differ significantly. Therefore, the average detection efficiency does not strongly depend on the K + K − invariant mass.
The detection efficiency is determined using both the φK + K − weighted PHSP MC and K + K − K + K − PHSP MC. The combined detection efficiency is given by where ǫ i and N i denote the detection efficiency and the signal yields of the ith mode, respectively. N total is the total signal yield obtained by fitting the K + K − K recoil-mass data, N 1 is the φK + K − signal yield, N 2 = N total − N 1 , and ǫ is the weighted detection efficiency for the K + K − K + K − final state. Figure 6 shows a comparison of the normalized momentum spectra of the kaon between the data and the weighted MC result for √ s = 3.080 GeV.

VI. DETERMINATION OF THE BORN CROSS SECTION
The Born cross section is calculated by where N obs is the number of observed signal events, L is the integrated luminosity, (1 + δ) stands for (1 + δ r ) · (1 + δ v ), and (1 + δ r ) is the ISR correction factor, which is obtained by a QED calculation [30] and by taking the line shape of the Born cross section measured by the BABAR experiment into account. The vacuum polarization factor (1+δ v ) is taken from a QED calculation with an accuracy of 0.5% [31], and ǫ is the detection efficiency. The branching fraction of the intermediate process φ → K + K − (49.2 ± 0.5%) [25] is taken into account in the determination of the cross section of e + e − → φK + K − .
Both ǫ and (1 + δ) are obtained from MC simulations of the signal reaction for each c.m. energy. In the conexc generator, the cross section for the ISR process (σ e + e − →γX ) is parametrized using where √ s ′ is the effective c.m. energy of the final state with s ′ = s(1 − x), x depends on the energy of the radiated photon according to x = 2E γ / √ s, W (s, x) is the radiator function and Π( √ s ′ ) describes the vacuum polarization (VP) effect. The latter includes contributions from leptons and quarks. The detection efficiency and the radiative-correction factor depend on the input cross section, and are determined by an iterative procedure, in which the line shape of the cross section from BABAR is used initially, and the updated Born cross section is obtained according to the simulation. We repeat the procedure until the measured Born cross section does not change by more than 0.5%.
The values of L, N obs , (1+δ) and ǫ are listed in Table I, together with the measured cross section at each energy point. Figures 7 (a) and 7 (b) show the line shapes of cross sections for e + e − → K + K − K + K − and e + e − → φK + K − , respectively.
The center-of-mass energy ( √ s), integrated luminosity (L), the yields of signal events (N obs ), the product of radiative correction factor and vacuum polarization factor (1 + δ), detection efficiency (ǫ), and Born cross section (σ B ). The first uncertainties are statistical and the second systematic.

VII. SYSTEMATIC UNCERTAINTY
Several sources of systematic uncertainties are considered in the measurement of the Born cross sections. These include the luminosity measurements, the differences between the data and the MC simulation for the tracking efficiency, PID efficiency, kinematic fit, the fit procedure, the MC simulation of the ISR-correction factor and the vacuum-polarization factor, as well as un-certainties in the branching fractions of the decays of intermediate states.
(a) Luminosity: The integrated luminosity of the data samples used in this analysis are measured using largeangle Bhabha scattering events, and the corresponding uncertainties are estimated to be 1.0% [18].
(b) Tracking efficiency: The uncertainty of the tracking efficiency is investigated using a control sample of the e + e − → K + K − π + π − process [32]. The difference in tracking efficiency between data and the MC simulation is estimated to be 1% per track. Hence, 3.0% is taken as the systematic uncertainty for the three selected kaons.
(c) PID efficiency: To estimate the uncertainty in the PID efficiency, we study K ± PID efficiencies with the same control samples as those used in the tracking efficiency. The average difference in PID efficiency between data and the MC simulation is found to be 1% per charged track. Therefore, 3.0% is taken as the systematic uncertainty for the three selected kaons.
(d) Kinematic fit: The uncertainty associated with the kinematic fits comes from the inconsistency of the track helix parameters between data and the MC simulation. The helix parameters for the charged tracks of MC samples are corrected to eliminate the inconsistency, as described in Ref. [33], and the agreement of χ 2 distributions between data and the MC simulation is significantly improved. We take the differences of the selection efficiencies with and without the correction as the systematic uncertainties.
(e) Fit procedure: A fit to mass spectrum of the recoiling kaon is performed to determine the signal yields of the e + e − → K + K − K + K − process, and the two kaon invariant mass M (K + K − ) is fitted to determine the number of e + e − → φK + K − events. The following three aspects are considered when evaluating the systematic uncertainty associated with the fit procedure.
(1) Fit range: The M (K ± ) spectrum of the recoiling kaon is fitted by varying the range from (0.3, 0.7) GeV/c 2 to (0.31, 0.69) GeV/c 2 . The M (K + K − ) spectrum is fitted in the region from 0.98 to 1.15 GeV/c 2 . An alternative fit range, from 0.98 to 1.20 GeV/c 2 , is considered. The differences between the yields are treated as the systematic uncertainty from the fit range.
(2) Signal shape: The signal shape of the mass spectrum of the recoiling kaon is described by a shape obtained from a MC simulation convolved with a Gaussian function. The uncertainty related to this line shape is estimated with an alternative fit using the same line-shape function, but fixing the width of the Gaussian function to a value differing by one standard deviation from the width obtained in the nominal fit. The signal shape of the φ is described by a P -wave Breit-Wigner function convolved with a Gaussian function. An alternative fit with a MC shape convolved with a Gaussian function is performed. The difference in yield between the various III: Relative systematic uncertainties (in %) for the cross section of e + e − → K + K − K + K − . The uncertainties are associated with the luminosity (L), tracking efficiency (Tracking), PID efficiency (PID), fit range (Range), signal and background shape (Sig. shape and Bck. shape), the initial-state radiation factor (ISR), the vacuum-polarization correction factor (VP), the weighted detection efficiency (ǫ), MC statistics (MC) and others. The total uncertainty is obtained by summing the individual contributions in quadrature.  Table III  fits is considered as the systematic uncertainty from the signal shape.
(3) Background shape: The background shape of the mass spectrum for the recoiling kaon is described as a second-order Chebyshev polynomial function. A fit with a first-order Chebyshev polynomial function for the background shape is used to estimate its uncertainty. The background shape for φ-mass distribution is described (f) ISR factor: The cross section is measured by iterating until (1+δ r )ǫ converges, and the difference between the last two iterations is taken as the systematic uncertainty associated with the ISR-correction factor.
(g) VP factor: The uncertainty on the calculation of the VP factor is 0.5% [31].
(h) Branching fraction: The experimental uncertainties in the branching fraction for the process φ → K + K − are taken from the PDG [25].
(i) Weighted detection efficiency: The detection efficiencies obtained in different processes are combined using the previously-described method. The combined uncertainty is calculated by accounting for the statistical variation, by one standard deviation, of the signal yields.
To obtain a reliable detection efficiency of e + e − → φK + K − , the PHSP MC sample is weighted to match the distribution of the background-subtracted data. To consider the effect on the statistical fluctuations of the signal yield in the data, a set of toy-MC samples, which are produced by sampling the signal yield and its statistical uncertainty of the data in each bin, are used to estimate the detection efficiencies.
(j) MC statistics: The uncertainty is estimated by the number of the generated events, whereby the weighting factor has been taken into account.
(k) Other systematic uncertainties: Other sources of systematic uncertainties include the trigger efficiency, the determination of the start time of an event, and the modeling of the final-state radiation in the simulation. The total systematic uncertainty due to these sources is estimated to be less than 1.0%. To be conservative, we take 1.0% as its systematic uncertainty.
Assuming all of the above systematic uncertainties, shown in Tables III and IV, are independent, the total systematic uncertainties are obtained by adding the individual uncertainties in quadrature.

VIII. SUMMARY AND DISCUSSION
In summary, using data collected with the BESIII detector taken at twenty c.m. energies from 2.100 to 3.080 GeV, we present measurements of the processes e + e − → K + K − K + K − and φK + K − and we obtain the corresponding Born cross sections. The Born cross sections of the process e + e − → K + K − K + K − are in good agreement with the results by BABAR, but with improved precision. The Born cross sections for the channel e + e − → φK + K − are measured for the first time at twenty energy points. Both data sets reveal anomalously high cross sections at √ s = 2.232 GeV.
A previous analysis on a much smaller dataset [17] has demonstrated that the K + K − K + K − final state exhibits resonant substructure. It is difficult to disentangle these contributions from other final states, and we make no attempt to do so.
By examining the φK + K − cross section as a function of c.m. energy, an enhancement at √ s = 2.232 GeV, i.e. near the ΛΛ production threshold, is observed. The cross section of e + e − → ΛΛ is also found to be anomalously high at the threshold [34]. In the case of charged baryons one would expect a Coulomb enhancement factor, which, however, is absent in the of the electricallyneutral Λ. It has been suggested that a narrow resonance, very close to the threshold, might provide an explanation [35]. BABAR has observed an enhancement at 2.175 GeV and a sharp peak at 2.3 GeV, corresponding to φK + K − final states with K + K − invariant masses smaller than 1.06 GeV/c 2 and within a mass interval of 1.06−1.2 GeV/c 2 , respectively. The intriguing φ(2170) resonance [13] has a relatively wide width and it is very close to the kinematical threshold, but not close enough to be related to the observed anomaly. Alternatively, the enhancement at 2.232 GeV could be explained by an interference effect of different resonances. More data in the vicinity would be helpful to understand the anomaly. (a) Comparison of the measured Born cross section of e + e − → K + K − K + K − to that of previous measurements [16]. The gray circles are from BABAR, the red rectangles are the results obtained in this work. The BESIII results include statistical and systematical uncertainties. The errors of the BABAR data only include the statistical uncertainty.
(b) Born cross section of e + e − → φK + K − obtained in this work. For BESIII data, the errors reflect both statistical and systematical uncertainties.