Observation of a resonant structure in $e^{+}e^{-} \to K^{+}K^{-}\pi^{0}\pi^{0}$

Using $e^{+}e^{-}$ collision data samples with center-of-mass energies ranging from 2.000 to 2.644 GeV, collected by the BESIII detector at the BEPCII collider, and with a total integrated luminosity of 300 pb^{-1}, a partial-wave analysis is performed for the process $e^{+}e^{-} \to K^{+}K^{-}\pi^{0}\pi^{0}$. The total Born cross sections for the process $e^{+}e^{-} \to K^{+}K^{-}\pi^{0}\pi^{0}$, as well as the Born cross sections f or the subprocesses $e^{+}e^{-} \to \phi \pi^{0}\pi^{0}$, $K^{+}(1460)K^{-}$, $K^{+}_{1}(1400)K^{-}$, $K^{+}_{1}(1270)K^{-}$ and $K^{*+}(892)K^{*-}(892)$, are measured versus the center-of-mass energy. The corresponding results for $e^{+}e^{-} \to K^{+}K^{-}\pi^{0}\pi^{0}$ and $\phi \pi^{0}\pi^{0}$ are consistent with those of BaBar and have much improved precision.By analyzing the cross sections for the four subprocesses, $K^{+}(1460)K^{-}$, $K^{+}_{1}(1400)K^{-}$, $K^{+ }_{1}(1270)K^{-}$ and $K^{*+}K^{*-}$, a structure with mass $M$ = (2126.5 $\pm$ 16.8 $\pm$ 12.4)~MeV/c^{2} and width $\Gamma$ = (106.9 $\pm$ 32.1 $\pm$ 28.1)~MeV is observed with an overall statistical significance of 6.3 $\sigma$, although with very limited significance in the subprocesses $e^{+}e^{-} \to K^{+}_{1}(1270)K^{-}$ and $K^{*+}(892)K^{*-}(892)$. The resonant parameters of the observed structure suggest it can be identified with the $\phi(2170)$, thus the results provide valuable input to the internal nature of the $\phi(2170)$.

Besides an excited φ state, the quark model also predicts excited ρ and ω states in the 2 GeV/c 2 mass range [23]. Finding this set of excited vector mesons would help establish the corresponding ρ, ω, and φ meson families and would set a baseline for theoretical models. Since these excited vector mesons can each decay into K ( * ) K ( * ) final states, analyzing the K ( * ) K ( * ) invariant mass spectra in e + e − annihilation becomes an effective means to discover them.
In this Letter, we present a PWA of the process e + e − → K + K − π 0 π 0 using data collected with the BE-SIII detector. The ten data samples used in this analysis have center-of-mass (c.m.) energies ranging from 2.000 to 2.644 GeV and have a total integrated luminosity of 300 pb −1 . The c.m. energy values and integrated luminosities of each data set are presented in Table I in the supplemental material [24]. Charge-conjugated processes are always included by default.
Detailed descriptions of the design and performance of the BESIII detector can be found in Ref. [25]. A Monte Carlo (MC) simulation based on Geant4 [26], including the geometric description of the BESIII detector and its response, is used to optimize the event selection criteria, estimate backgrounds, and determine the detection efficiency. The signal MC samples are generated using the package ConExc [27], which incorporates a higherorder ISR correction. Background samples of the processes e + e − → e + e − , µ + µ − and γγ are generated with the Babayaga [28] generator, while e + e − → hadrons and two photon events are generated by the Luarlw [29] and Bestwogam [30] generators, respectively.
The selection criteria for charged tracks, particle identification (PID), and photon candidates are the same as those in Ref. [31].
The process e + e − → K + K − π 0 π 0 results in the final state K + K − γγγγ. Thus, candidate events with only two oppositely-charged kaons and at least four photons are selected. To improve the kinematic resolution and suppress background, a six-constraint (6C) kinematic fit imposing energy-momentum conservation, as well as two additional π 0 mass constraints, is carried out under the hypothesis e + e − → K + K − π 0 π 0 . The combination with minimum χ 2 6C is retained for further analysis. The candidate events are required to satisfy χ 2 6C < 80. After the above selection criteria, detailed studies indicate that the backgrounds are negligible.
Using the GPUPWA framework [32], a PWA is performed on the surviving candidate events to disentangle the intermediate processes present in e + e − → K + K − π 0 π 0 . The quasi two-body decay amplitudes in the sequential decays are constructed using covariant tensor amplitudes [33]. The intermediate states are parameterized with relativistic Breit-Wigner (BW) functions, except for the f 0 (980), which is described with a Flatté formula [34]. The resonance parameters of the f 0 (980) and the wide resonance σ in the fit are fixed to those in Ref. [34] and Ref. [34,35], respectively, and those of other intermediate states are fixed to PDG values, or measured in the analysis. To include the resolution for the narrow φ(1020) resonance, a Gaussian function is convolved with the BW function, but this is not done for the other resonances. The relative magnitudes and phases of the individual intermediate processes are determined by performing an unbinned maximum likelihood fit using MI-NUIT [36].
We start the fit procedure by including all possible intermediate states in the PDG that conserve J PC , where these intermediate states can decay into K + K − , π 0 π 0 , K ± π 0 , K + K − π 0 or K ± π 0 π 0 final states. Then we examine the statistical significance of the individual amplitudes, and drop the ones with statistical significance less than 5σ. The process is repeated until no amplitude remains with a statistical significance less than 5σ. After that, all the removed processes are reintroduced individually to make sure that they are not needed in the fit. In the above approach, the statistical significance of each individual amplitude is determined by the changes in the negative log likelihood (NLL) value and the number of free parameters in the fit with and without the corresponding amplitude included.
The above strategy is performed individually on the data sets at √ s = 2.125 and 2.396 GeV, which have the largest luminosities among the ten data sets.
The nominal solution for data at √ s = 2.125 GeV includes the two-body decay 1270) and ω(1420)π 0 , as well as the three-body decay processes K + K − σ, For the data at √ s = 2.396 GeV, the additional intermediate processes K * + 2 (1430)K * − (892), K * + (892)K − π 0 and φ(1020)f 0 (1370) are included, but without the φ(1020)σ and φ(1020)f 2 (1270) processes. An interesting decay mode K * + (1410)K − , which is expected to have a sizeable decay rate for a conventional 3 3 S 1 ss state [9], is found to be less than 3σ in both data samples. In the above, the three-body decays are treated as consecutive quasi twobody decays with a very broad resonance decaying into K + K − or K + π 0 and modeled as a 1 − phase space distribution. The intermediate states K + (1460), K + 1 (1270), K + 1 (1400) decay into K * + (892)π 0 , and ω(1420) decays into K * ± (892)K ∓ , followed by K * + (892) → K + π 0 . The state K * + 0 (1430) decays into K + π 0 . The state φ(1020) decays into K + K − and σ, f 0 (980), f 2 (1270), f 0 (1370) decay into π 0 π 0 . The masses and widths of the K(1460), K 1 (1400), K 1 (1270) and ω(1420) in the fit are determined by scanning the likelihood value, and the results are consistent with the parameters in the PDG. The masses and widths of other intermediate states are fixed to PDG values. The statistical significance of all intermediate processes are summarized in sections II and III of the supplemental material [24], respectively. The corresponding comparison of invariant mass spectra and angular distributions between data and MC projections are shown in section IV of the supplemental material.
For the other eight data samples, due to limited statistics, we do not perform the above optimization strategy to determine which intermediate processes to include. Instead, we use the same intermediate processes as the data sets with nearby c.m. energy. The data sets with √ s = 2.000, 2.100, 2.175, 2.200 and 2.232 GeV (referred to as group I data), use the same processes as √ s = 2.125 GeV, while the other three points (group II data) use the same processes as √ s = 2.396 GeV. The total Born cross sections for e + e − → K + K − π 0 π 0 and the Born cross sections for the intermediate processes are obtained at each c.m. energy using: where N sig is the corresponding signal yield, and is determined by calculating the fraction according to the PWA results for the individual intermediate process; L int is the integrated luminosity; (1 + δ) r is the ISR correction factor obtained from a QED calculation [27,37] and incorporating the input cross section from this analysis iteratively; 1 |1−Π| 2 is the vacuum polarization factor taken from a QED calculation [38]; ǫ is the detection efficiency obtained from a PWA-weighted MC sample; and Br is the product of branching ratios of the intermediate states as quoted in the PDG [1]. In the decay e + e − → K + (1460)K − , the branching fraction of K(1460) → K * (892)π is included in the measured cross section since it has never been measured.
Two categories of systematic uncertainties are considered in the measurement of the Born cross sections. The first category includes uncertainties associated with the luminosity, track detection, PID, kinematic fit, ISR correction, and the branching fractions of intermediate states. The uncertainty associated with the integrated luminosity is 1% at each energy point [39]. The uncertainty of the detection efficiency is 1% for each charged track [40] and photon [41]. The PID efficiency uncertainty is 1.0% for each charged track [40]. The uncertainty related to the kinematic fit is estimated by correcting the helix parameters of the simulated charged tracks to match the resolution [42]. The uncertainty associated with the ISR correction factor is estimated to be the difference of (1 + δ r )ǫ between the last two iterations in the cross section measurement. The systematic uncertainties from the branching ratios of intermediate states in the subsequent decays are taken from the PDG [1]. The second category of uncertainties are from the PWA fit procedure. Fits with alternative scenarios are performed, and the changes of signal yields are taken as systematic uncertainties. Uncertainties from the BW parameterization are estimated by replacing the constant-width BW with the mass-dependent width. Uncertainties associated with the resonance parameters, which are taken from the PDG and fixed in the fit, are estimated by alternative fits superposing additional constraints on these resonance parameters, where the superposed constraints follow Gaussian distributions with widths equal to their uncertainties. One thousand fits are performed, and the resultant standard deviations of the signal yields are taken as systematic uncertainties. Uncertainties associated with the additional resonances are estimated by alternative fits including the components K * (1410)K or the K * 2 (1430)K * , which are most significant, but less than 5σ. Uncertainties due to the barrier factor are estimated by varying the radius of the centrifugal barrier from 0.7 to 1.0 fm. To estimate the uncertainties on the detection efficiency related to the fit parameters in the PWA, one hundred MC samples are generated with five hundred groups of parameters of PWA amplitudes which is sampled from a multi-variable Gaussian function according to their mean values and their covariance error matrix from the nominal fit. The standard deviations of the resultant detection efficiencies are considered as the uncertainties.
In the above procedure, the uncertainties associated with the barrier factor, resonance parameterization and additional resonances are strongly affected by the statistics. Thus, those uncertainties of data with √ s=2.125 GeV are assigned to the group I data, while those of data with √ s=2.396 GeV are assigned to the group II data. Assuming all sources of systematic uncertainties are independent, the total uncertainties are the quadratic sums of the individual values, shown in section V of the supplemental material [24], where the sources of the uncertainties tagged with '*' are assumed to be 100% correlated among each energy points.
The measured total Born cross sections for e + e − → K + K − π 0 π 0 and the Born cross sections for the subprocess e + e − → φπ 0 π 0 , summing over all the π 0 π 0 intermediate processes and their interferences, are shown in Fig. 1. Good agreement is found with the previous results from BaBar. In order to study the properties of 1 −− states, the cross sections for the processes e + e − → K + (1460)K − , K + 1 (1400)K − , K + 1 (1270)K − and K * + (892)K * − (892), referred to as the KK processes, are shown in Fig. 2. A clear peak between 2.1 and 2.2 GeV is present in the process e + e − → K + (1460)K − , and dips are observed for the processes e + e − → K + 1 (1400)K − and K + 1 (1270)K − in almost the same energy region. This may be due to destructive interference between different components. No obvious structure or dip is present in the process e + e − → K * + (892)K * − (892). All the various numbers used in the cross section calculation are summarized in section I of the supplemental material [24]. To further examine the structure, a binned χ 2 fit, incorporating the correlated and uncorrelated uncertainties among different energy points, is performed to the cross sections for the K + (1460)K − , K + 1 (1400)K − , K + 1 (1270)K − and K * + (892)K * − (892) processes. The fit probability density function (PDF) for the individual processes is the coherent sum of a continuum component f 1 and a resonant component f 2 : where φ is the relative phase between the two components. By considering phase space Φ( √ s), the energydependent cross section of the QED process, and the relative orbital angular momentum L in the two-body decay, the amplitude f 1 is described as where q is the momentum of the daughter particle. The resonant amplitude f 2 is described with a BW function, where M R is the mass of the structure, Γ R is the total width, Γ e + e − R is its partial width to e + e − , Br is the decay branching fraction to a given final state, and q 0 is the momenta of the daughter particle in the rest frame of the parent particle (M R ).
A simultaneous fit, assuming the same structure among the K + (1460)K − , K + 1 (1400)K − , K + 1 (1270)K − and K * + (892)K * − (892) processes, is performed to the measured cross sections, as shown in Fig. 2. In the fit, M R and Γ R are shared parameters between the four processes and are floated, while n, the production BrΓ e + e − R , and the relative phase angle φ are floated and final state dependent. For e + e − → K + 1 (1270)K − and K + 1 (1400)K − , L = 0, while L = 1 for the other two modes. The fit results have two solutions with equal fit quality, identical M R = (2126.5 ± 16.8) MeV/c 2 and Γ R = (106.9±32.1) MeV, but different BrΓ e + e − R and φ for the processes e + e − → K + 1 (1400)K − and K + 1 (1270)K − , as summarized in Table I. The statistical significance of the structure is estimated with the change of χ 2 (∆χ 2 ) and the number of degrees of freedom (∆ndof) between the scenarios with and without the structure included in the fit. The overall statistical significance is 6.3σ, obtained with ∆χ 2 =63.8 and ∆ndof=10. The significance of the resonant state for each KK process is also estimated and summarized in Table I. The significances of the resonant state in the processes e + e − → K + (1460)K − and K + 1 (1400)K − are greater than 4.5σ, while no significant signal is found in the other two processes. We also estimate the upper limit at the 90% confidence level on the production BrΓ e + e − R to be 1.9 eV for e + e − → K * + (892)K * − (892) and 12.5(297.6) eV for e + e − → K + 1 (1270)K − . The systematic uncertainties on the resonant parameters come from the absolute c.m. energy measurement, the measured cross section, and the fit procedure. The uncertainty of the c.m. energy from BEPCII is small, and is ignored in the determination of the parameters of the structure. The statistical and systematic uncertainties of the measured cross section are incorporated in the fit, thus no further uncertainty is necessary. The uncertainties associated with the fit procedure include those from the fit range and signal model. The uncertainty from the fit range is investigated by excluding the last energy point √ s = 2.644 GeV in the fit. The resultant changes, 5.1 MeV/c 2 for mass and 9.1 MeV for width, are taken as the systematic uncertainties. To assess the systematic uncertainty associated with the signal model, an alternative BW function with energy-dependent width is implemented in the fit, and results in differences of 11.3 MeV/c 2 and 26.5 MeV for mass and width, respectively, which are taken as the systematic uncertainties. The overall systematic uncertainties are the quadratic sum of the individual ones, 12.4 MeV/c 2 and 28.1 MeV for the mass and width, respectively.
In summary, a PWA of the process e + e − → K + K − π 0 π 0 is performed for ten data samples with c.m. energies from 2.000 to 2.644 GeV and with an integrated luminosity of 300 pb −1 . The Born cross sections for e + e − → K + K − π 0 π 0 and φπ 0 π 0 are obtained and are consistent with those from the BaBar experiment. We also measure the cross sections for the processes e + e − → K + (1460)K − , K + 1 (1400)K − , K + 1 (1270)K − , and K * + (892)K * − (892), individually, and perform a simultaneous fit on the obtained results. The fit results in a structure with mass M = (2126.5 ± 16.8 ± 12.4) MeV/c 2 , width Γ = (106.9 ± 32.1 ± 28.1) MeV, and statistical significance 6.3 σ, where the uncertainties are statistical and systematic, respectively. The structure is directly produced in e + e − collisions, thus has J P C = 1 −− . This structure has a mass close to the masses of the vector particles φ(2170), ρ(2150) and ω(2290) listed in the PDG [1]. Its width is only consistent with the φ(21770) and is different from the others by more than 3σ.