Measurement of cross sections of the interactions $e^+e^-\rightarrow \phi\phi\omega$ and $e^+e^-\rightarrow \phi\phi\phi$ at center-of-mass energies from 4.008 to 4.600 GeV

Using data samples collected with the BESIII detector at the BEPCII collider at six center-of-mass energies between 4.008 and 4.600 GeV, we observe the processes $e^+e^-\rightarrow \phi\phi\omega$ and $e^+e^-\rightarrow \phi\phi\phi$. The Born cross sections are measured and the ratio of the cross sections $\sigma(e^+e^-\rightarrow \phi\phi\omega)/\sigma(e^+e^-\rightarrow \phi\phi\phi)$ is estimated to be $1.75\pm0.22\pm0.19$ averaged over six energy points, where the first uncertainty is statistical and the second is systematic. The results represent first measurements of these interactions.


Introduction
The experimental understanding of hadron production in electron-positron annihilation has been achieved with the measurement of the total inclusive hadronic cross sections, the so-called R measurement [1], and the exclusive measurement of final states involving pions, kaons and other light hadrons at various centerof-mass (c.m.) energies [2,3]. The tools for describing the e + e − annihilation to hadrons process generally include the use of the KKMC generator [4], which includes initial and final state radiation, and the Pythia [5] program based on the Lund String model or Parton Shower model that hadronizes the final-state quarks. The KKMC-Pythia combination is not expected to correctly describe the processes with more than two vector mesons in the final state, as they correspond to higher order Quantum Chromodynamics (QCD) processes and are generally associated with multiple gluons. The experimental results provide more constraints Preprint submitted to Physics Letters B March 19, 2018 on the higher-order QCD calculation.
The BaBar and Belle collaborations reported the observation of significant double charmonium production e + e − → J/ψcc and found the ratio σ(e + e − → J/ψcc)/σ(e + e − → J/ψX) to be ∼ 0.6 [6], which indicates that a surprisingly large fraction of e + e − → J/ψX events are produced by the e + e − → J/ψcc process. This experimental result has stimulated much theoretical interest. Various theoretical approaches, such as NRQCD factorization [7] and the light cone method [8], have been proposed to make corrections to the low ratio predicted by the non-relativistic calculation, which predicts a much lower value for the cross section [9]. The validity of the theoretical investigations can be tested over a wide kinematical range with double or triple quarkonia (ss, cc, bb) produced in e + e − annihilations. In particular, strangeonia ss are located in the region of transition between perturbative QCD and nonperturbative QCD. The e + e − annihilation to multiple ss states may provide an important experimental opportunity in the low-energy region.
In this paper, we report on the first measurement of the Born cross sections of e + e − → φφω and e + e − → φφφ processes at c.m. energies E cm = 4.008, 4.226, 4.258, 4.358, 4.416 and 4.600 GeV [10]. The data samples were collected by the BESIII detector at the BEPCII collider [11].
Additionally, we also measure the ratio σ(e + e − → φφω)/σ(e + e − → φφφ), where many of the systematic uncertainties are canceled. The mixing angle of ω and φ is expected to be small and its effect on the ratio can be neglected. In the e + e − annihilation process, without considering the intermediate resonance, the final φφφ states would be generated via one virtual photon and two gluons or three virtual photons, as illustrated in FIG. 1. The production via two virtual photons and one gluon is forbidden, because the gluon carries color while the final state is color neutral. By replacing s(7)s(8) with (uū + dd)/ √ 2 in Fig. 1(a), we obtain the ratio σ(e + e − →γ * gg→2(ss)+(uū+dd)/ √ 2) σ(e + e − →γ * gg→3(ss)) ∼ ( 4 9 + 1 9 )/2 1 9 = 2.5, because the vertex "A" is proportional to the charge squared of the quarks. If, on the other hand, (uū + dd)/ √ 2 is substituted for s(3)s(4) or s(5)s(6), the ratio would be about 1 since the strong interaction vertex only relies on the mass of the quarks. Considering the above two cases in Fig. 1(a) and neglecting the small contribution from Fig. 1(b), σ(e + e − →γ * gg→2(ss)+(uū+dd)/ √ 2) σ(e + e − →γ * gg→3(ss)) would range from 1 to 2.5, depending on the ratio of the two cases above. The study of σ(e + e − → φφω)/σ(e + e − → φφφ) can therefore help to understand the production mechanism of e + e − annihilation to multiple quarkonia.

Detector and Monte Carlo Simulation
The BESIII detector, as described in detail in Ref. [12], has a geometrical acceptance of 93% of the solid angle. A small-cell, helium-based main drift chamber (MDC) immersed in a 1 T magnetic field measures the momentum of charged particles with a resolution of 0.5% at 1 GeV/c, and provides energy loss (dE/dx) measurements with a resolution better than 6% for electrons from Bhabha scattering. The electromagnetic calorimeter (EMC) detects photons with a resolution of 2.5% (5%) at an energy of 1 GeV in the barrel (end cap) region. A time-of-flight system (TOF) assists in particle identification (PID) with a time resolution of 80 ps (110 ps) in the barrel (end cap) region.
A GEANT4-based [13] Monte Carlo (MC) simulation software, which includes the geometric description of the BESIII detector and the detector response, is used to optimize the event selection criteria, determine the detection efficiency and estimate background contributions. The simulation includes the beam energy spread and initial-state radiation (ISR) modeled with KKMC [4]. In this analysis, 0.5 million events of e + e − → φφω and e + e − → φφφ are generated individually at different c.m. energies corresponding to the experimental values. Both processes are simulated with a uniform distribution in phase space (PHSP). The observed cross sections for e + e − → φφω and e + e − → φφφ at the six energy values in this analysis are used as the inputs in the KKMC simulation for ISR effects. In line with the partial reconstruction technique that is implemented in the analysis, the signal process e + e − → φφω is simulated with both φ decaying into K + K − and the ω decaying into all possible final states, while in the simulation of e + e − → φφφ events, all three φ are generated to decay via all possible modes.

Event Selection
The candidate events for e + e − → φφω and φφφ are selected with a partial reconstruction method to get higher efficiencies. We reconstruct two φ mesons with their prominent K + K − decay mode and identify the remaining ω or φ meson with the mass recoiling against the reconstructed φφ system.
For each charged track, the polar angle in the MDC must satisfy | cos θ| < 0.93, and the point of closest approach to the e + e − interaction point must be within ±10 cm in the beam direction and within 1 cm in the plane perpendicular to the beam direction. We identify charged kaon candidates using the dE/dx and TOF information. The probabilities L(π) and L(K) are determined for the π and K hypothesis, respectively. Kaons are identified by requiring L(K) > L(π).
The φ candidates are formed from pairs of identified kaons with opposite charges. Their invariant mass is required to satisfy 1.01 < M (K + K − ) < 1.03 GeV/c 2 . At least two φ candidates with no shared tracks are required in each event. If there are more than two φ candidates in one event, only the φφ combination with the minimum ∆M is kept for further analysis, and the two φ candidates are randomly labeled as φ 1 or φ 2 . The mass difference ∆M is defined as is the nominal mass of the φ meson taken from the particle data group (PDG) [14]. Fig. 2 (a) depicts the scatter plot of M φ1 (K + K − ) versus M φ2 (K + K − ) by combining the data samples at six c.m. energies. A clear accumulation of events is observed around the intersection of the φ 1 and φ 2 mass regions, which indicates e + e − → φφX signals. The mass of the system recoiling against the reconstructed φφ is where E cm is the c.m. energy obtained by analyzing the di-muon process e + e − → γ ISR/FSR µ + µ − , with a precision of 0.02% [10]. E φφ and p φφ are the energy and momentum of the reconstructed φφ pair in the e + e − rest system. As shown by the solid points in Fig. 2 (b), we obtain two clear peaks in the vicinities of ω and φ in the RM (φφ) distribution, which indicates the processes e + e − → φφω and φφφ, respectively.

Study of Backgrounds in RM (φφ)
To ensure that the observed ω and φ signal in the RM (φφ) distribution originate from the processes e + e − → φφω and φφφ, we perform a study of the potential peaking backgrounds. The two dimensional (2D) sidebands illustrated in Fig. 2 (a) are used to study the potential background without a φφ pair in the final state, where the φ sidebands are defined as 0.99 < M (K + K − ) < 1.00 GeV/c 2 and 1.04 < M (K + K − ) < 1.06 GeV/c 2 . The non-φ 1 and/or nonφ 2 processes are estimated by the weighted sum of the events in the horizontal and vertical sideband regions, with the entries in the diagonal sidebands subtracted to compensate for the double counting of the background without any φ in final state. The weighting factor for the φ 2 but non-φ 1 events in the horizontal sidebands is the ratio of the number of φ 2 but non-φ 1 events under the signal region (n sig bkg ) to the number of φ 2 but non-φ 1 events in the horizontal sidebands (n sdb bkg ). n sig bkg and n sdb bkg are determined from the 2D fit to The weighting factor for the φ 1 but nonφ 2 (non-φ 1 and non-φ 2 ) events in the vertical (diagonal) sidebands are determined similarly. The 2D probability density functions for the components φ 1 φ 2 , φ 1 but nonφ 2 , non-φ 1 but φ 2 , non-φ 1 and non-φ 2 are constructed by the product of two one-dimensional functions. The φ peak is described with a MC-derived shape convoluted with a Gaussian function to take into account the resolution difference between data and MC simulation. The non-φ component is described with second-order polynomial functions. The estimated RM (φφ) distribution with weighted 2D sidebands events is shown as  the open circles in Fig. 2 (b). Since the φ signal is close to the K + K − production threshold, we are not able to obtain a sideband which is far enough away from the signal region at the lower side of M (K + K − ). Thus, the small ω and φ signals observed in RM (φφ) estimated with the 2D sideband are from the leakage of the real e + e − → φφ + ω/φ signals. From studies of signal MC samples, the ratio of the signal events in the 2D sideband regions to those in the signal region is estimated to be 3%∼5%. We also estimate the peaking background in the RM (φφ) distribution for the process e + e − → φφφ with the MC samples. The dominant peaking backgrounds is from the e + e − → K + K − φφ and e + e − → K + K − K + K − φ processes. When the directly produced K + K − (K + K − K + K − ) is reconstructed as φ (φφ), these two processes would contribute as peaking backgrounds in the RM (φφ) distribution. The contamination rate of the e + e − → K + K − φφ (e + e − → K + K − K + K − φ) events to e + e − → φφφ is estimated to be ∼ 1.0% (0.1%) at each energy point with the assumption that the c.m. energy dependent cross section for e + e − → K + K − φφ (e + e − → K + K − K + K − φ) is the same as for e + e − → φφφ.
We take 1.0% as the uncertainty on the size of the peaking backgrounds of e + e − → φφφ. Similarly, the dominant peaking backgrounds of e + e − → φφω is from the e + e − → K + K − φω and e + e − → K + K − K + K − ω processes. For e + e − → φφω, the uncertainty from the peaking backgrounds is determined to be 1.0%.

Fits to the RM (φφ) Spectrum and Cross Section Results
The reconstruction efficiencies and yields of e + e − → φφω and φφφ signals are determined by the fit to the RM (φφ) distribution for MC simulation and data, respectively.

Correction to RM (φφ)
Compared with the values in the PDG, the measured masses of the ω and φ mesons in the RM (φφ) distribution deviate to the left with ∼4.5 MeV. This deviation may be induced by ISR, the energy loss of the reconstructed kaons and final state radiation (FSR), or the uncertainty of E cm . The overall effect is considered as a shift on E cm , ∆E cm .
We estimate ∆E cm by studying the process e + e − → φK + K − with partially reconstructing one φ meson and one charged kaon. The recoil mass against the reconstructed φK is calculated with RM (φK) = (E cm − E φK ) 2 − p 2 φK , where E φK and p φK are the energy and momentum of the reconstructed φK in the system of e + e − . ∆E cm is estimated with is approximately m(K) from PDG and E φK is the average over all φK + K − events. RM (φφ) for each event is then corrected by subtracting ∆RM (φφ) in the data and MC samples, where ∆RM (φφ) = Ecm−E φφ RM(φφ) × ∆E cm . As a consequence, the measured masses of the ω and φ mesons obtained by fitting the RM (φφ) distributions are consistent with the values in the PDG.

Fits to the RM (φφ) Spectrum
An unbinned maximum likelihood fit is performed to the corrected RM (φφ) distributions. The signal distribution is modeled by the MC-derived signal shape. The study of the selected φ signal indicates that the mass resolution difference for the φ signal is very small. Therefore, we assume the resolution of RM (φφ) is the same between data and MC simulation, and the corresponding systematic uncertainty will be considered. The background shape is described by a third-order Chebyshev polynomial function with parameters fixed to the values obtained by fitting all samples together, since some samples have small statistics. The corresponding fit results are shown in Fig. 3. The statistical significances of the ω/φ signals are examined using the differences in likelihood values of fits with and without an ω/φ signal component included in the fits. Both ω and φ signals are seen with statistical significances of more than 3σ for each data sample, and the significances of ω and φ are both larger than 10 σ if all six data samples are combined. The yields of ω and φ signal events and the corresponding statistical significances for each sample are summarized in Table 1 and Table 2, respectively.

Reconstruction Efficiency
The e + e − → φφω and φφφ signal MC samples are simulated by assuming a uniform distribution in phase space. The reconstruction efficiency of the two reconstructed φs depends on their production angles. The comparison of the cosine of the polar angles θ for the two reconstructed φ mesons between data and MC simulation is presented in Fig. 4, where the cos θ distributions are obtained by fitting the RM (φφ) distribution for events with cos θ in given bins. All the data samples are combined, assuming the cos θ distributions do not depend on the c.m. energy. To take into account the deviation in cos θ distributions between the data and the PHSP MC samples, the reconstruction efficiencies are determined with PHSP MC samples incorporating the re-weighting correction according to the 2D distribution of cos θ 1 versus cos θ 2 of data and PHSP MC samples.

Cross Section Results
The Born cross section is calculated by where N obs is the number of observed signal events, L int is the integrated luminosity, (1+δ r ) is the radiative correction factor, (1 + δ v ) is the vacuum polarization factor, ǫ is the detection efficiency including reconstruction and all selection criteria, and B is the branching fraction of φ → K + K − . The vacuum polarization factor is taken from a QED calculation. With the input of the observed c.m. energy dependent σ(e + e − → φφω) and σ(e + e − → φφφ), and using a linear interpolation to obtain the cross sections in the full range, the radiative correction factor is calculated in QED [15].
Since the radiative correction factor and the detection efficiency both depend on the line shape of the input cross section, the Born cross sections of e + e − → φφω and e + e − → φφφ are determined with four iterations until convergence has been reached. The values of all variables used in the calculation of σ(e + e − → φφω) and σ(e + e − → φφφ) are listed in Table 1 and Table 2, respectively.

Systematic Uncertainties of Cross Sections
Several sources of systematic uncertainties are considered in the measurement of the Born cross sections. These include differences between the data and the MC simulation for the tracking efficiency, PID efficiency, mass window requirement, the MC simulation of the radiative correction factor and the vacuum polarization factor. We also consider the uncertainties from the fit procedure, the peaking backgrounds, the simulation model as well as uncertainties of the branching fraction of φ → K + K − and the integrated luminosity.
a. Tracking efficiency. The difference in tracking efficiency for the kaon reconstruction between the data and the MC simulation is estimated to be 1.0% per  Summary of the measurements of the e + e − → φφω process. Listed in the table are the c.m. energy Ecm, the integrated luminosity L int , the number of the observed events N obs , the reconstruction efficiency ǫ, the vacuum polarization factor (1 + δ v ), the radiative correction factor (1 + δ r ), the measured Born cross section σ B , and statistical significance. The first uncertainty of the Born cross section is statistical, and the second is systematic.  Summary of the measurements of the e + e − → φφφ process. Listed in the table are the c.m. energy Ecm, the number of the observed events N obs , the reconstruction efficiency ǫ, the radiative correction factor (1 + δ r ), the measured Born cross section σ B , and statistical significance. The first uncertainty of the Born cross section is statistical, and the second is systematic. The integrated luminosity L int and the vacuum polarization factor (1 + δ v ) are same with those in Table 1.   track [16]. Therefore, 4.0% is taken as the systematic uncertainty for four kaons.
b. PID efficiency. PID is required for the kaons, and the uncertainty is estimated to be 1.0% per kaon [16]. Hence, 4.0% is taken as the systematic uncertainty of the PID efficiency for four kaons.
c. φ mass window. A mass window requirement on the K + K − invariant mass might introduce a systematic uncertainty on the efficiency. The reconstructed φ signals are fit with a MC shape convoluted with a Gaussian function that describes the disagreement between data and MC simulation. The mean and width of the Gaussian function are left free in the fit, which turn out to be close to 0 within 3 times of uncertainty. The systematic uncertainty from the M (K + K − ) requirement is ignored.
d. Fit procedure. For the six data samples, the yields of e + e − → φφω and φφφ events are obtained by a fit to the distribution of the mass recoiling against the reconstructed φφ system. The following two aspects are considered when evaluating the systematic uncertainty associated with the fit procedure. (1) Signal shape.-In the nominal fit, the signal shapes are described by the MC shape obtained from MC simulation. An alternative fit with the MC shape convoluted with a Gaussian function for the ω/φ signal shape is performed, where the parameters of the Gaussian function are free. The resulting difference in the yield with respect to the nominal fit is considered as the systematic uncertainty from the signal shape. This uncertainty is negligible compared to the statistical uncertainty.
(2) Background shape.-In the nominal fit, the background shape is described with a third-order Chebyshev polynomial function. The fit with a fourth-order Chebyshev polynomial function for the background shape is performed to estimate the uncertainty due to the background parametrization.
e. Peaking backgrounds. The uncertainty is taken as 1.0%, as described in Sec. 4.
f. Line shape of cross section. The line shape of the e + e − → φφω and φφφ cross sections affects the radiative correction factor and the reconstruction efficiency. The corresponding uncertainty is estimated by changing the input of the observed line shape within one standard deviation.
g. vacuum polarization factor. The QED calculation used to determine the vacuum polarization factor has an accuracy of 0.5% [17].
h. Simulation model. The differences between the efficiencies obtained with and without re-weighting the PHSP MC sample are taken as the uncertainties associated with the simulation model.
i. Luminosity. The time-integrated luminosity [18] of each sample is measured with a precision of 1% with Bhabha events.
j. Branching fractions. The uncertainty in the branching fraction for the process φ → K + K − is taken from the PDG [14].
Assuming all of the systematic uncertainties shown in Tables 3 and 4 are independent, the total systematic uncertainties are obtained by adding the individual uncertainties in quadrature.
7. Ratio σ(e + e − → φφω)/σ(e + e − → φφφ) The right plot of Fig. 5 shows the measured ratios r cs ≡ σ(e + e − → φφω)/σ(e + e − → φφφ) at different c.m. energy, and the statistical-weighted average. Except for the measurement at 4.008 GeV, the ratios are consistent with each other within one statistical standard deviation. In the calculation of r cs , many uncertainties on the cross sections cancel, such as the uncertainties in the tracking, PID, B(φ → K + K − ) and luminosity. Only the uncertainties from the background shape, line shape and MC simulation model are considered in the determination of r cs . From the measurements at six energy points in Table 5, we obtain the statistical-weighted average r cs = 1.75 ± 0.22 ± 0.19, where the first uncertainty is statistical and the second systematic. The systematic uncertainties of r cs at different c.m. energies are assumed to be independent in this calculation.

Summary and Discussion
With the data samples collected between 4.008 and 4.600 GeV with the BESIII detector, the processes e + e − → φφω and e + e − → φφφ are observed for the first time. The Born cross sections are determined at six c.m. energies and the average ratio σ(e + e − → φφω)/σ(e + e − → φφφ) over the six c.m. energies is calculated to be 1.75 ± 0.22 ± 0.19, which is in the range of the estimation with Fig. 1. Our measurements of these two processes provide experimental constraints on the theoretical calculations of the three vectors production in the e + e − annihilation.

Acknowledgement
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National Key Basic Research Program of China under Contract No. 2015CB856700; National Natural Science Foundation of China (NSFC) under Contracts Nos. 11235011,