Cross sections for the reactions $e^+ e^-\to K_S^0 K_L^0$, $K_S^0 K_L^0 \pi^+\pi^-$, $K_S^0 K_S^0 \pi^+\pi^-$, and $K_S^0 K_S^0 K^+K^-$ from events with initial-state radiation

We study the processes $e^+ e^-\to K_S^0 K_L^0 \gamma$, $K_S^0 K_L^0 \pi^+\pi^-\gamma$, $K_S^0 K_S^0 \pi^+\pi^-\gamma$, and $K_S^0 K_S^0 K^+K^-\gamma$, where the photon is radiated from the initial state, providing cross section measurements for the hadronic states over a continuum of center-of-mass energies. The results are based on 469 fb$^{-1}$ of data collected with the BaBar detector at SLAC. We observe the $\phi(1020)$ resonance in the $K_S^0 K_L^0$ final state and measure the product of its electronic width and branching fraction with about 3% uncertainty. We present a measurement of the $e^+ e^-\to K_S^0 K_L^0 $ cross section in the energy range from 1.06 to 2.2 GeV and observe the production of a resonance at 1.67 GeV. We present the first measurements of the $e^+ e^-\to K_S^0 K_L^0 \pi^+\pi^-$, $K_S^0 K_S^0 \pi^+\pi^-$, and $K_S^0 K_S^0 K^+K^-$ cross sections, and study the intermediate resonance structures. We obtain the first observations of \jpsi decay to the $K_S^0 K_L^0 \pi^+\pi^-$, $K_S^0 K_S^0 \pi^+\pi^-$, and $K_S^0 K_S^0 K^+K^-$ final states.

We study the processes e + e − → K 0 S K 0 L γ, K 0 S K 0 L π + π − γ, K 0 S K 0 S π + π − γ, and K 0 S K 0 S K + K − γ, where the photon is radiated from the initial state, providing cross section measurements for the hadronic states over a continuum of center-of-mass energies.The results are based on 469 fb −1 of data collected with the BABAR detector at SLAC.We observe the φ(1020) resonance in the K 0 S K 0 L final state and measure the product of its electronic width and branching fraction with about 3% uncertainty.We present a measurement of the e + e − → K 0 S K 0 L cross section in the energy range from 1.06 to 2.2 GeV and observe the production of a resonance at 1.67 GeV.We present the first measurements of the e + e − → K 0 S K 0 L π + π − , K 0 S K 0 S π + π − , and K 0 S K 0 S K + K − cross sections, and study the intermediate resonance structures.We obtain the first observations of J/ψ decay to the K 0 S K 0 L π + π − , K 0 S K 0 S π + π − , and K 0 S K 0 S K + K − final states.

I. INTRODUCTION
The idea to use electron-positron annihilation events with initial-state radiation (ISR) to study processes with energies below the nominal e + e − center-of-mass (E c.m. ) energy was outlined in Ref. [1].The possibility of exploiting ISR to measure low-energy cross sections at high-luminosity φ and B factories is discussed in Refs.[2][3][4], and motivates the study described in this paper.This is of particular interest because of a threestandard-deviation discrepancy between the current measured value of the muon anomalous magnetic moment (g − 2) and that predicted by the Standard Model [5], in which hadronic loop contributions are obtained from experimental e + e − annihilation cross sections at low E c.m. energies.The study of ISR events at B factories provides independent results over a continuum of energy values for hadronic cross sections in this energy region and also contributes to the investigation of low-mass resonance spectroscopy.
Studies of the ISR processes e + e − → µ + µ − γ [6, 7] and e + e − → X h γ, using data of the BABAR experiment at SLAC, where X h represents any of several exclusive multihadron final states, have been reported previously.The studied final states include: charged hadron pairs π + π − [7], K + K − [8], and pp [9]; four or six charged mesons [10][11][12]; charged mesons plus one or two π 0 mesons [11][12][13][14]; and a K 0 S plus charged mesons [15].Together, these demonstrate good detector efficiency for events of this kind, and well understood tracking, particle identification, and π 0 and K 0 S reconstruction.This paper reports measurements of the K 0 S K 0 L , K 0 S K 0 L π + π − , K 0 S K 0 S π + π − , and K 0 S K 0 S K + K − final states, produced in conjunction with a hard photon, that is assumed to result from ISR. Candidate K 0 S decays are reconstructed in the π + π − decay mode.This is the first ISR measurement from BABAR that includes K 0 L mesons, which we detect via their nuclear interactions in the electromagnetic calorimeter.We use the e + e − → γφ → γK 0 S K 0 L reaction to measure the K 0 L detection efficiency directly from the data.The e + e − → K 0 S K 0 L cross section is measured from threshold to 2.2 GeV.For the other final states, we measure cross sections from threshold to 4 GeV, study the internal structure of the events, and perform the first measurements of their J/ψ branching fractions.Together with our previous measurements [8,11], these results provide a much more com-plete understanding of the KK, KKππ, and KKKK final states in e + e − annihilations.

II. THE BABAR DETECTOR AND DATA SET
The data used in this analysis were collected with the BABAR detector at the PEP-II asymmetric-energy e + e − storage ring.The total integrated luminosity used is 468.6 fb −1 [16], which includes data collected at the Υ (4S) resonance (424.7 fb −1 ) and at a c.m. energy 40 MeV below this resonance (43.9 fb −1 ).
The BABAR detector is described in detail elsewhere [17].Charged particles are reconstructed using the BABAR tracking system, which comprises the silicon vertex tracker (SVT) and the drift chamber (DCH) inside the 1.5 T solenoid.Separation of pions and kaons is accomplished by means of the detector of internally reflected Cherenkov light and energy-loss measurements in the SVT and DCH.The hard ISR photon, photons from π 0 decays, and K 0 L are detected in the electromagnetic calorimeter (EMC).Muon identification, provided by the instrumented flux return, is used to select the µ + µ − γ final state.
In order to study the detector acceptance and efficiency, we have developed a special package of simulation programs for radiative processes based on the approach suggested by Kühn and Czyż [18].Multiple collinear soft-photon emission from the initial e + e − state is implemented with the structure-function technique [19,20], while additional photon radiation from the final-state particles is simulated using the PHOTOS package [21].The precision of the radiative simulation does not contribute more than 1% uncertainty to the efficiency calculation.
The four-meson final states are generated according to a phase-space distribution.We simulate the K 0 S K 0 L γ channel using a model that includes the φ(1020) and two additional resonances, fitted to all available e + e − → K 0 S K 0 L cross section measurements [22][23][24][25][26], which cover the range from threshold up to about 2.5 GeV.Samples of roughly five times the number of expected events are generated for each final state and processed through the detector response simulation [27].These events are then reconstructed using the same software chain as the data.Variations in detector and background conditions are taken into account.
We also simulate a number of background processes.Based on our experience with final states including kaons, we consider the ISR processes K 0 S K 0 L (φ)η, K 0 S K ± π ∓ , K 0 S K ± π ∓ π 0 , K 0 S K 0 L π + π − π 0 , and K 0 S K 0 L π 0 , with normalizations based on our previous measurements and isospin relations.In addition, we generate a large sample of the as-yet unobserved final state K 0 S K 0 L π 0 π 0 γ, which is a potential background.We also simulate several non-ISR backgrounds, including e + e − → qq (q = u, d, s, c) events using the Jetset7.4 [28] generator, and e + e − → τ + τ − events using the KORALB [29] generator.

III. THE ISR PHOTON AND K 0
S SELECTION Photons are reconstructed as clusters of energy deposits in contiguous crystals of the EMC.We consider the cluster in the event with the highest energy in the e + e − c.m. frame, and require ISR event candidates to contain a cluster with E γ c.m. > 3 GeV, which we denote as the ISR photon.The ISR photon detection efficiency has been studied using µµγ events [7], and we apply a polar-angle-dependent correction of typically −1.5±0.5% to the simulated efficiency.
In these events, we reconstruct K 0 S candidates decaying to two charged pions from pairs of oppositely charged tracks not identified as electrons.They must have a well reconstructed vertex between 0.2 and 40.0 cm in radial distance from the beam axis, and their total momentum must be consistent with originating from the interaction region.The m(π + π − ) invariant mass distribution for these K 0 S candidates is shown in Fig. 1 for both data (points) and a simulation (histogram) containing only genuine K 0 S .The background level is relatively low and we select candidates in the 482 < m(π + π − ) < 512 MeV/c 2 mass range (vertical lines on Fig. 1), and use the sidebands 472-482 and 512-522 MeV/c 2 to estimate the contribution from non-K 0 S backgrounds.A few thousand events (about 1% of the total number of events) have more than one selected K 0 S candidate, and we use these to study the K 0 S K 0 S π + π − and K 0 S K 0 S K + K − final states.Considering only the "best" K 0 S candidate, with m(π + π − ) closest to the nominal [30] K 0 S mass, we also include these events in the K 0 S K 0 L and K 0 S K 0 L π + π − measurements.The K 0 S detection efficiency has been studied very carefully at BABAR, with data-MC differences in the efficiency determined as a function of the K 0 S direction and momentum.We apply a correction event by event, which introduces an overall correction +1.1±1.0% to the number of K 0 S .We also require the event to contain exactly zero or two tracks that are consistent with originating from the interaction region, excluding those in the selected K 0 S candidate(s).Any number of additional tracks and EMC clusters is allowed.
The decay length of the K 0 L meson is large, and the probability to detect a K 0 L decay in the DCH is low.Instead, we look for a cluster in the EMC resulting from the interaction of a K 0 L with a nucleus in the EMC material.Such clusters are indistinguishable from photon-induced clusters, and give poor resolution on the K 0 L energy.In this section, we describe the use of a clean sample of e + e − → φγ, φ → K 0 S K 0 L events to optimize our selection of K 0 L clusters and measure their detection efficiency and angular resolution.In sections VI and VII, we describe the use of the selected K 0 L candidate clusters to study the φ resonance and measure the e + e − → K 0 S K 0 L cross section above the φ region, respectively.
Using the four-momenta of the best selected K 0 S , the ISR photon, and the initial electron and positron, we can calculate the recoil-mass-squared, where E 0 = E + + E − and p 0 = p + + p − are the energy and total momentum vector of the initial e + e − system, E γ and p γ (with E γ ≡ | p γ |) are the energy and momentum vector of the photon, and E K 0 S and p K 0 S are the energy and momentum vector of the K 0 S candidate.The presence of the reaction e + e − → K 0 S K 0 L γ would be evident as a peak in the m rec distribution at the mass of the K 0 L .Due to the large uncertainty of the measured ISR photon energy, the calculated value of m rec also has a large uncertainty.However, if we assume the reaction e + e − → γφ(1020) → γK 0 S K 0 L , we can calculate the constrained ISR photon energy E c γ according to: where n γ is a unit vector along the ISR photon direction and m φ is the φ meson mass [30].Using E c γ instead of the measured E γ in Eq. 1, we obtain a much better resolution on the recoil mass m c rec for genuine events of that type.The m c rec distribution for our data is shown in Fig. 2 as the points.A simulated distribution for genuine e + e − → γφ → γK 0 S K 0 L events is shown as the histogram.Selecting events with m c rec > 0.4 GeV/c 2 (corresponding to m(K 0 S K 0 L ) < 1.1 GeV/c 2 ), with the additional requirement that there be no other track within a 0.2 cm radius of the interaction point, we obtain a very clean sample of K 0 S K 0 L γ events, without any need to detect the K 0 L meson.
The non-K 0 S background, estimated from the sidebands of the m(π + π − ) distribution in Fig. 1, contributes 0.8% of the events in Fig. 2.This background arises from e + e − → γγ events in which one photon converts to a misidentified electron-positron pair.We estimate backgrounds from other ISR final states containing a real K 0 S using the simulation.Normalized contributions to the m c rec distribution for 3, cumulatively, as shaded, hatched, and open histograms.The simulated backgrounds from e + e − → qq (q = u, d, s, c), and e + e − → τ τ events are found to be negligible.Fitting the simulated non-K 0 S and ISR backgrounds with smooth functions and summing them together with the signal simulation, we obtain excellent agreement with the observed spectrum.The total background is 6.9±0.5% of the selected events.
The position of the K 0 L peak in Fig. 3 is very sensitive to both the reconstructed K 0 S candidate mass and the assumed φ-meson mass (see Eqs. 1,2).There is a small 0.21±0.02MeV/c 2 data-MC difference in the K 0 S peak position in Fig. 1.As a cross check, we correct the data for this difference and vary m φ in Eq. 2 for the data so that the experimental m c rec peak position matches The simulated signal distribution is shown as the dashed histogram and the sum of all simulated events as the solid histogram.
that of the simulation.This results in an estimate of m φ = 1019.480±0.040±0.036MeV/c 2 , where the systematic uncertainty includes the effects of the nominal K 0 mass (0.024 MeV/c 2 [30]), the K 0 S momentum measurement in the DCH (0.020 MeV/c 2 ), and the DCH-EMC mis-alignment (0.018 MeV/c 2 ).This is consistent with the value, tabulated by the Particle Data Group (PDG) m φ = 1019.455± 0.020 MeV/c 2 [30].
Subtracting the non-K 0 S and ISR-produced backgrounds, we obtain 81 012 ± 285 (447 434 for the MC simulation) K 0 S K 0 L γ events in the φ mass region without requiring K 0 L detection.These events must satisfy our trigger and software filters, which were designed for various classes of events.We study efficiencies in data and simulation using prescaled events not subject to these filters, and obtain a correction of (+3.9 ± 2.3)%.Furthermore, the pions from K 0 S decays in this particular reaction have a relatively large probability to overlap in the DCH, and the reconstruction efficiency for overlapping tracks is not well simulated.We introduce a +1.5±0.6% correction for this effect.

B. The K 0 L detection efficiency
We select events with m c rec > 0.47 MeV/c 2 (vertical line in Fig. 2), reducing the background level from 6.9% to 2.8%.Using the K 0 S and ISR photon angles and momenta, we calculate the hypothetical K 0 L direction for each event, and look for an EMC cluster in that direction.∆ψ, rad.

Energy of cluster, GeV
FIG. 4: The EMC cluster energy versus the opening angle between the measured cluster direction and the predicted K 0 L direction for all non-ISR clusters in the data.Events/0.02rad.
FIG. 5: The difference in azimuthal angle between the EMC cluster direction and the predicted K 0 L direction for selected clusters in the data (points) and MC-simulated Figure 4 shows a two-dimensional plot of the EMC cluster energy versus the opening angle ∆ψ between the predicted K 0 L direction and measured cluster direction for all clusters in the data except those assigned to the ISR photon.A clean signal is observed at high cluster energies, but the background from low-energy clusters is large.We consider clusters with energy greater than 0.2 GeV, and select the one closest to the predicted K 0 L direction if it is within 0.5 radians.This yields K 0 L detection probabilities of about 48% in the data and 51% in the MC simulation.
We then study the resolution in polar (θ) and azimuthal (φ) angle of the selected K 0 L clusters as a function of their position in the detector and the predicted K 0 L energy.The resolutions in the two angles are consistent, with no significant dependence on position or energy.The overall ∆φ distributions are shown for data and simulation in Fig. 5. Good agreement is seen, with root-mean-square deviations of 0.035 radians.We use this value in the kinematic fits.

V. THE KINEMATIC FIT PROCEDURE
Each candidate event selected in Sec.III is subjected to a set of constrained kinematic fits in which the fourmomenta and covariance matrices of the initial e + e − , the ISR photon, the best K 0 S candidate, and the two tracks from the interaction region, if present, are taken into account.The three-momentum vectors for each particle including the photon obtained from these fits are determined with better accuracy, and are used in further calculations.
First, we consider each neutral cluster with E > 0.2 GeV (excluding the ISR photon) as a K 0 L candidate, and perform a three-constraint (3C) kinematic fit under the The angular resolutions for K 0 L clusters discussed in the previous section are used, and the K 0 L momentum is determined in the fit.We retain the K 0 L candidate cluster giving the best χ 2 value in each event.
We then perform a kinematic fit under the K 0 S K ± π ∓ π 0 γ hypothesis, where the cluster is assumed to be one photon from a π 0 decay, rather than a K 0 L .Such events can enter the sample if a charged kaon is misidentified as a pion and only one photon from the π 0 decay is considered.Similarly, we perform fits under the hypotheses of the other backgrounds discussed in Sec.II, giving us additional χ 2 variables with which to suppress these processes.
We perform additional fits to the events with more than one K 0 S candidate under the For each pair of K 0 S candidates, a four-constraint (4C) kinematic fit is performed using the four-momenta and covariance matrices of all initial-and final-state particles.The combination with the best χ 2 for each hypothesis is retained.

A. Additional selection criteria and background subtraction
To study this mass region, we consider events selected as described in Sec.IV A, with m c rec > 0.4 GeV/c 2 (see Fig. 2).We select a K 0 L cluster where possible, using the 3C fits described in Sec.V, and obtain the χ 2 distribution for the best K 0 S K 0 L γ candidate shown in Fig. 6 as the points.The unshaded histogram is for the corresponding MC-simulated pure K 0 S K 0 L γ events, normalized to the data in the region χ 2 < 10 where we expect very low background.FIG.6: The three-constraint χ 2 distribution for K 0 S -cluster-γ events in the data (points) fitted under the K 0 S K 0 L γ hypothesis.The open histogram represents the same distribution for the MC-simulated signal events, normalized to the data in the region χ 2 < 10, and the shaded histogram represents the estimated background in the data.
The experimental and simulated distributions are broader than a typical 3C χ 2 distribution due to multiple soft-photon emission from the initial state, which is not taken into account in the fit, but is present in both the data and simulation.The observed difference at higher χ 2 values is due to background in the data and possibly a data-MC difference in the angular uncertainty of the K 0 L cluster.
For further analysis we require χ 2 (K 0 S K 0 L ) < 15 (vertical line in Fig. 6), and for these events we calculate the K 0 L candidate mass according to Eqs. 1,2 and perform the background subtraction described in Sec.IV A. We obtain 27 925 ± 176 events for the data (871 background events are subtracted) and 164 179 events for the MCsimulation, representing samples with the K 0 L detected.Dividing by the corresponding numbers of events before K 0 L cluster selection, we obtain K 0 L detection efficiencies, including the effects of the kinematic fit and χ 2 selection, of 0.3447±0.0017for the data and 0.3724±0.0008for the simulation.The double ratio 0.9394±0.0052 is applied as a correction factor to account for this data-MC difference.This ratio is independent of momentum and polar angle of the K 0 L .We use the four-vectors returned by the kinematic fit to calculate the K 0 S K 0 L invariant mass, the distribution of which is shown in Fig. 7.The φ(1020) resonance is clearly visible, with a width of about 10 MeV, much larger than the nominal width of the resonance [30] due to the resolution of this final state.The background, estimated as described above, is shown as the shaded histogram.We fit it with a smooth, empirical function, shown as the line, and use the fit result in each bin for background subtraction.

B. Fit for the φ(1020) parameters
To obtain the parameters of the φ(1020), we fit the background-subtracted distribution in Fig. 7 with a cross section σ(s) convolved with a resolution matrix Res(j, i).
In each mass bin i: where efficiency correction factor for the χ 2 cut, track overlap, and event filter, L(s) is the ISR luminosity, calculated at leading-order [4], and N 0 (j) is the acceptance-corrected number of events expected for bin j.
The 100×100 resolution matrix is obtained from simulation by binning the reconstructed versus simulated K 0 S K 0 L invariant mass for signal events in 1×1 MeV/c 2 intervals.The distribution of differences between the reconstructed and simulated masses near 1.020 GeV/c 2 , corresponding to a row of this matrix, is shown in Fig. 8.The K 0 S K 0 L threshold and a radiative tail are visible.We normalize each row to unit area, and introduce an additional variable Gaussian smearing σ add , to account for any data-MC difference in the resolution.
We describe the cross section near the φ resonance using formulae discussed in detail in Refs.[22,31]: where is the propagator for a vector resonance V , C = 0.389 × 10 12 nb MeV 2 /c 4 [30], and describes the energy-dependent width, and for the φ we use the set of final states f and ηγ, with corresponding branching fractions B(V → f ) and phase space factors P V →f (s).We include the influence of the ρ(770) and ω(782) resonances in the in the energy-dependent width according to the "ideal" quark model, which assumes their decay rates to K 0 S K 0 L are a factor of two lower than that of the φ.We use the relation in Eq. 4 for the corresponding cross sections.We introduce a complex constant A K 0 S K 0 L to describe the contributions of higher radial excitations of the ρ, ω and φ mesons to the cross section, as well as any deviations from the "ideal" quark structure relations for the ρ(770) and ω(782).It can be written in terms of two free parameters, a non-resonant cross section σ bkg , and a phase Ψ, The fitted value of Ψ is consistent with zero, and we fix it to zero in the final fit, but propagate its fitted uncertainty as a systematic uncertainty to account for model dependence.The detection efficiency, shown as a function of mass in Fig. 9, is obtained by dividing the number of selected MC-simulated events in each 0.001 GeV/c 2 mass interval by the number generated in the same interval.The mass dependence is well described by a linear fit, which we use in all calculations.This efficiency includes the geometrical acceptance of the detector for the final-state photon and the charged pions from the K 0 S decay, the inefficiency of the detector subsystems, and event losses due to additional soft-photon emission from the initial state.It is not sensitive to the detector mass resolution.
The result of the fit is projected on the backgroundsubtracted invariant mass distribution in Fig. 10.We obtain the resonance parameters: where the first uncertainties are statistical, the second systematic, and the third due to model dependence, evaluated by varying σ bkg by its uncertainty.
We introduce an additional Gaussian smearing to describe an uncertainty in the detector resolution, and obtain σ add = 0.6 ± 0.2 MeV/c 2 , which improves the χ 2 of the fit in the 1.0-1.05GeV/c 2 region from 59 to 53, for 51 degrees of freedom.We estimate systematic uncertainties of 0.05 MeV/c 2 in mass and 0.05 MeV in width from the uncertainty of the σ add value.The other systematic uncertainties are summarized in Table I, along with the corrections applied to the measurements.A total correction of +14.1±2.9% is applied to the number of events.The largest contribution to the uncertainty is from the software filter, due to the limited number of available prescaled events.
Our parameter values are consistent with the most precise cross section measurement, σ φ = 1376 ± 24 nb [22], and with the PDG values m φ = 1019.455± 0.020 MeV/c 2 and Γ φ = 4.26 ± 0.04 MeV [30].Since each row of the resolution matrix is normalized to unit area, the smearing procedure does not affect the total number of events, which is proportional to the product Γ φ σ 0 of the total width and peak cross section of the φ.Using this product as a free parameter in the fit, we obtain the product of the electronic width of the φ and its branching fraction to 0019 keV, where the first uncertainty is statistical, the second systematic, and the third due to model dependence.Using B K 0 S K 0 L = 0.342 ± 0.004 or Γ φ = 4.26 ± 0.04 MeV from Ref. [30], we obtain Γ ee = 1.228 ± 0.037 ± 0.014 keV or B ee B K 0 S K 0 L = 0.986 ± 0.030 ± 0.009, respectively, where the first uncertainty is our total experimental uncertainty and the second is from the PDG tables.These values are consistent with the most recent measurement of Γ ee = 1.235±0.022keV [32], and with the PDG values of Γ ee = 1.27±0.04keV and B ee B K 0 S K 0 L = 1.006±0.016[30], and have comparable precision.
In this section we consider events with m(K 0 S K 0 L ) > 1.06 GeV/c 2 .Since the e + e − → K 0 S K 0 L cross section drops much more rapidly with increasing mass than the background, we apply additional selection criteria compared to the criteria of Sec.VI A. In all cases, we consider K 0 S candidates with 0.482 < m(π + π − ) < 0.512 MeV/c 2 (see Fig. 1) and use sideband data to subtract non-K 0 S background from all studied quantities.

A. Additional selection criteria
We consider all EMC clusters except those assigned to the ISR photon and the K 0 L as photon candidates, and combine each pair into a π 0 candidate.Figure 11 shows a scatter plot of the higher of the energies E γ max of the two photons assigned to the pair versus the corresponding diphoton mass m γγ .A large signal from events containing a π 0 is observed.To reduce this background, we require E γ max < 0.5 GeV (horizontal line in Fig. 11).
Since signal event may contain several background clusters, this reduces the signal efficiency.We measure this loss using events in the φ region, where no π 0 signal is observed, but background clusters are present in both data and simulation.We find losses of 10% in the data and 7% in the simulation, and apply the 3% difference as a correction.The 3C χ 2 distribution for the remaining candidate events with m(K 0 S K 0 L ) > 1.06 GeV/c 2 is shown as the points in Fig. 12.The open histogram shows the corresponding simulated distribution for genuine K 0 S K 0 L γ events, normalized to the data in the region χ 2 < 3. The shaded, cross-hatched and hatched areas represent the simulated contributions from the ISR channels φη, K 0 S K 0 L π 0 , and K 0 S K 0 L π 0 π 0 , respectively.These channels contribute significant background, and almost entirely account for the difference between the data and signal-MC χ 2 distributions.We find no significant contribution from simulated non-ISR backgrounds.
We select events with χ 2 (K 0 S K 0 L ) < 10, and use events from the control region 10 < χ 2 (K 0 S K 0 L ) < 20 (vertical lines in Fig. 12) to estimate the background in the signal region.The signal region contains 6264 data and 13 292 MC-simulated events, while the control region contains 2968 and 2670, respectively.

B. Background subtraction
To obtain any distribution of the K 0 S K 0 L signal events N d 0 (m), we take the experimental events in the signal Events/unit χ 2

FIG. 12:
The 3C χ 2 distributions for K 0 S K 0 L γ candidate events in the data (points) and signal simulation (open histogram), fitted under the K 0 S K 0 L hypothesis.The shaded, cross-hatched, and hatched areas represent the simulated contributions from the ISR φη, K 0 S K 0 L π 0 , and K 0 S K 0 L π 0 π 0 channels, respectively.region of Fig. 12, N d s (m), and subtract the background events, taken from the control region N d c (m), corrected for the presence of signal events, estimated from MCsimulation N MC c (m): where b = 1.15 is the simulated ratio of background events in the signal and control regions, and a = N d 0 (m)/N MC s (m) is a factor equalizing the number of signal and simulated events.This procedure relies on good agreement between data and simulation in both the χ 2 and mass distributions.As noted above, the MC simulation uses a "world average" cross section, well measured below 1.4 GeV, but based on only the DM1 measurement [24], which has large statistical uncertainties, in the 1.4-2.4GeV E c.m. region.We adopt an iterative procedure, in which we reweight the simulated mass distribution to match our measurement and repeat the subtraction until there is no change in the results.
Figure 13a shows the K 0 S K 0 L mass distribution for data events in the χ 2 control region of Fig. 12 as points, with the shaded histogram showing the distribution for signal MC at the final iteration.The signal contribution is not large, and the difference between the data and weighted MC-simulated distributions is scaled by b = 1.15 to estimate the background in the χ 2 signal region.The squares in Fig. 13b represent this background estimate in each K 0 S K 0 L invariant mass bin.We also estimate the background directly from the MC simulation of the ISR φη, K 0 S K 0 L π 0 , and K 0 S K 0 L 2π 0 processes, shown as the histogram in Fig. 13b.The two estimates agree relatively well, but the MC simulation does not incorporate the correct mass distributions for these processes, and other unknown processes might contribute.The mass distribution after background subtraction is shown in Fig. 13c.In the 1.4-2.4GeV/c 2 mass region we select about 1000 events, compared with only 58 events found by the DM1 [24] experiment.

C. Simulated detection efficiency
The selection procedures applied to the data are also applied to the MC-simulated event sample.The resulting K 0 S K 0 L invariant mass distribution is shown in Fig. 14(a) for the signal and control (shaded histogram) regions.The mass dependence of the detection efficiency is obtained by dividing the number of reconstructed MC events in each mass interval by the number generated in that interval.The results are shown in Fig. 14(b).The 40 MeV/c 2 mass intervals used are wider than the detector resolution of 10 MeV/c 2 , but a small effect of the resolution on the efficiency is visible, due to the very steep decrease in the cross section with increasing mass.
D. The e + e − → K 0 S K 0 L cross section for c.m. energies above 1.06 GeV The cross section for e + e − annihilation into K 0 S K 0 L can be calculated from ) is the corresponding detection efficiency, estimated from the MC simulation with corrections to the K 0 S and ISR photon detection efficiencies; and the factor ǫ corr takes into account the data-MC differences in K 0 L detection and background filter efficiencies, and the energy requirement on additional photon clusters.The radiative correction R is within one percent of unity, with an estimated precision of about 1%.The differential luminosity dL(E c.m. ) associated with the interval dE c.m. centered at an effective collision energy of E c.m. is calculated using the leading order formula (see for example Ref. [13]), and the systematic uncertainty associated with the luminosity determination is estimated to be 0.5%.
The cross section for the reaction e + e − → K 0 S K 0 L after all corrections is shown as a function of energy in Fig. 15 and listed in Table II.Only statistical uncertainties are shown.The systematic uncertainty is dominated by the background subtraction procedure and is strongly correlated across the entire energy range.It is about 10% at 1.6 GeV, and increases with decreasing cross section to  12).The histogram in (a) represents MC-simulated signal events in the control region, and that in (b) represents the simulated background from φη, K 0 S K 0 L π 0 , and K 0 S K 0 L 2π 0 events; The squares show the total estimated background, obtained as the difference between the data and simulated distributions in (a), normalized to the data in the signal region.(c) The K 0 S K 0 L invariant mass distribution above 1.06 GeV/c 2 for the data after background subtraction.[23] and SND [25] experiments at VEPP-2M, and is much more precise than the result from the OLYA experiment [26].
In the 1.4-2.4GeV region, our result is much more precise than the only other available measurement, from the DM1 [24] experiment.
The measured cross section exhibits a distinctive structure around 1.6 GeV, indicating the presence of a vector resonance, perhaps the φ(1680).Denoting it φ ′ , we fit the cross section using Eq. 4 with the additional amplitude The energy-dependent width (see Eq. 5) assumes the branching fractions and phase space factors of the major φ(1680) decay modes, f = K * K, φη, φππ, and K 0 S K 0 L taken from Refs.[11,15].We fix the φ(1020) parameters to the values obtained in Sec.VI B, and float the parameters of the φ ′ .Since the other vector states in this energy range, such as ω(1420, 1650) and ρ(1450, 1700), are relatively wide and overlap considerably, we again describe the sum of their contributions using the nonresonant cross section σ bkg and phase Ψ of Eq. 7. First, we fix both to zero, and the fit yields a relatively good description of the data, with χ 2 =30 for (29-4) degrees of freedom.The result of the fit (solid curve) is compared with the data in Figs.16 and 17.
Next, we allow σ bkg and Ψ to float in the fit, and obtain Ψ = 0.2 ± 0.6 radians.Since this is consistent with zero, we fix it to zero and repeat the fit.The result is shown as the dashed curves in Figs.16 and 17.We obtain an improved description of the cross section, with χ 2 = 21 for (29-6) degrees of freedom and the fitted parameter values: where the first uncertainties are statistical and the second systematic, dominated by the difference between fixed and floated Ψ.The relative phase between the nonresonant background and the φ resonance is consistent with that between the φ ′ and φ resonances, but the un- certainty is very large.
Our parameter values for this resonance are consistent with those of the PDG for the φ(1680), and with the results of similar fits performed in Refs.[11,15] for the K * K, φη, and φππ decay modes of the φ(1680).However, as shown in Fig. 17, the cross section for e + e − → K + K − is quite different from that for K 0 S K 0 L , indicating substantial interference between the iso-scalar and iso-vector amplitudes in this energy range.The fitting function used above is not able to reproduce the FIG. 16: The e + e − → K 0 S K 0 L cross section (points) compared with the results of the fits described in the text with the non-resonant amplitude fixed to zero (solid lines) and floating (dashed lines).K + K − data [8], and therefore the results should be taken with caution.A simultaneous fit to the cross sections for e + e − → π + π − (pure iso-vector), π + π − π 0 (pure isoscalar), K + K − , K 0 S K 0 L , and perhaps other multihadron final states is needed to extract the iso-scalar and isovector components, together with reliable resonance parameter values.
The product Γ φ ′ σ φ ′ is proportional to the total number of events and does not depend on the experimental resolution.Using this product as a free parameter in the fit and Eq. 6 we obtain for the φ(1680) candidate where the first uncertainty is statistical, the second systematic, and the third due to model dependence.There is no independent measurement of the φ(1680) → K 0 S K 0 L branching fraction that could be used to calculate Γ ee .However, we have also measured Γ ee B K * K = 369±53 eV, Γ ee B φη = 138±43 eV [15], and Γ ee B φππ = 42±5 eV [11].
We assume these are the four dominant decay modes, estimate their rates, and use them in our Γ(s) calculation.
We now consider the events with exactly two tracks not assigned to the K 0 S candidate, but consistent with originating from the same event vertex.This final state has four charged particles, and therefore large backgrounds from ISR and non-ISR multihadron events.We make additional requirements on the two tracks and the rest of the event in order to suppress these backgrounds.

A. Additional selection criteria
The two additional tracks must not be identified as K ± , and are required to extrapolate to within ±3 cm of the collision point in the direction along the beam axis and 0.25 cm in the perpendicular direction.The event must contain no other tracks that extrapolate to within 1 cm of the axis, which is also the lower limit on the radial position of the K 0 S → π + π − vertex.Considering all pairs of EMC clusters except those assigned to the ISR photon and K 0 L candidates, we observe a large signal from π 0 , similar to that shown in Fig. 11.As in that case, we require E γ max < 0.5 GeV, reducing backgrounds from several sources with a loss of 3% in signal efficiency, as shown in Sec.VII A.
ISR K 0 S K ± π ∓ π 0 events with the charged kaon misidentified as a pion and a cluster from a π 0 photon taken as the K 0 L candidate are indistinguishable from signal events.To reduce this background, we pair the K 0 L cluster with all other EMC clusters.For every such pair with M (γγ) within 0.03 GeV/c 2 of the π 0 mass, we perform a kinematic fit to the K 0 S K ± π ∓ π 0 γ hypothesis and require χ 2 (K 0 S Kππ 0 ) > 100.The 3C χ 2 distribution for the remaining candidate events under the K 0 S K 0 L π + π − γ hypothesis is shown as the points in Fig. 18a, with the corresponding MC-simulated pure K 0 S K 0 L π + π − γ events shown as the open histogram.The simulated distribution is normalized to the data in the region χ 2 < 1, where the contribution of higher-order ISR is small and the background contamination is lowest, but still amounts to about 15% of the signal.The shaded, cross-hatched, and hatched areas represent the simulated contributions from the ISR φη, ISR K 0 S K ± π ∓ π 0 , and non-ISR q q channels, respectively.These backgrounds account for only half of the observed data-MC difference in the distribution at large χ 2 values.
We define a signal region χ 2 (K 0 S K 0 L π + π − ) < 25 and a control region 25 < χ 2 (K 0 S K 0 L π + π − ) < 50 (vertical lines in Fig. 18), from which we estimate backgrounds in the signal region.The signal region contains 10 788 data and 6825 MC events, while the control region contains 5756 and 633 events, respectively.

B. Background subtraction
The background to the K 0 S K 0 L π + π − mass distribution is subtracted in two stages.The χ 2 distributions for the K 0 S K ± π ∓ π 0 and non-ISR q q events peak at low values, since their kinematics are similar to those of signal events.We therefore subtract their MC-simulated contribution from both the signal and control regions of Fig. 18(a).There are large uncertainties in their normalizations, but this has little effect on the total uncertainty.The mass distribution for the data in the signal region before background subtraction is shown in Fig. 18(b) as the points, with the simulated K 0 S K ± π ∓ π 0 and q q events shown as the shaded and cross-hatched histograms, respectively.
We estimate the remaining background using the mass distributions for the remaining events in the signal and control regions, according to Eq. 8 of Sec.VII B. The contribution is shown as the hatched area in Fig. 18(b).We fit the sum of all backgrounds with a polynomial function to reduce the statistical fluctuations (curve in Fig. 18(b)) and use this fit for the background subtraction.The resulting mass distribution for e + e − → K 0 S K 0 L π + π − events is shown in Fig. 18(c).We observe 3320 events in the mass range from threshold to 4.0 GeV/c 2 .In addition to a main peak around 2 GeV/c 2 , a J/ψ signal and a possible structure just below 3 GeV/c 2 are visible.
We estimate the systematic uncertainty due to the background subtraction to be about 10% for m(K 0 S K 0 L π + π − ) < 2.5 GeV/c 2 (i.e., a 30% uncertainty on a 30% total background), increasing to about 30% in the 2.5-3.0GeV/c 2 region and reaching 100% above 3.4 GeV/c 2 , where background dominates.

C. Simulated detection efficiency
The selection procedures applied to the data are also applied to the MC-simulated event sample.The resulting K 0 S K 0 L π + π − invariant mass distribution is shown in Fig. 19(a) for the signal and control (shaded histogram) regions.The detection efficiency as a function of mass is obtained by dividing the number of reconstructed MC events in each 0.05 GeV/c 2 mass interval by the number generated in that interval, and is shown in Fig. 19(b).The 50 MeV/c 2 mass interval used is wider than the detector resolution of about 25 MeV/c 2 .Since the cross Events/0.05GeV/c 2 FIG.18: (a) The three-constraint χ 2 distributions for data (points) and MC simulated K 0 S K 0 L π + π − γ events (open histogram).The shaded, cross-hatched, and hatched histograms represent the simulated contributions from ISR φη, ISR K 0 S Kππ 0 , and non-ISR q q events, respectively.(b) The K 0 S K 0 L π + π − invariant mass distribution for data events in the signal region of (a) (points).The shaded and cross-hatched histograms represent the simulated contributions from ISR φη+K 0 S Kππ 0 and non-ISR q q events, respectively, and the hatched area represents that estimated from the control region.The curve shows the empirical fit used for background subtraction.(c) The K 0 S K 0 L π + π − invariant mass distribution after background subtraction.section has no sharp structures (except for the J/ψ signal, which is discussed below), we apply no corrections for the resolution.We apply all the corrections discussed above for data-MC differences in the tracking, photon, and K 0 L detection efficiencies.
FIG. 20: The e + e − → K 0 S K 0 L π + π − cross section.The error bars are statistical only.

D. The e
The cross section for the reaction e + e − → K 0 S K 0 L π + π − is calculated using Eq. 9 with the corrections described above, plus an additional 3% correction for the requirement on the maximum energy of extra EMC clusters.The cross section is shown as a function of energy in Fig. 20, and listed in Table III.There are no previous measurements for this final state.The cross section shows a threshold rise at 1.5 GeV, a maximum value of about 1 nb near 2 GeV, and a slow decrease toward higher energies, perturbed by the J/ψ signal.Only statistical uncertainties are shown.The total systematic uncertainty is dominated by the background subtraction procedure.It amounts to about 10% at 2 GeV, where the cross section peaks, and increases with decreasing cross section to ∼30% near 1.5 and 3 GeV, and to 100% well above 3 GeV.We fit these projections with a sum of two Breit-Wigner functions and a function describing the nonresonant contribution, yielding 3335 ± 115 K * (892) ± → K 0 S π ± decays, 3200 ± 151 K * (892) ± → K 0 L π ± decays, and a total of 286 ± 99 K * 2 (1430) ± decays.The total number of K * (892) ± decays is larger than the number of K 0 S K 0 L π + π − events, indicating correlated production of K * (892) + K * (892) − pairs.In each 0.04 GeV/c 2 bin of K 0 L π ∓ mass we fit the K 0 S π ± mass distribution with the same function, and the resulting numbers of K * (892) ± decays are shown in Fig. 23.
A strong signal of 2098 ± 61 ± 200 K * (892) ± is observed, where the second uncertainty is due to variations of the fitting procedure.This corresponds to the production of K * (892) + K * (892) − pairs in about 63% of all observed K 0 S K 0 L π + π − events.We also find 105 ± 23 ± 50 events at the K * 2 (1430) ± mass, corresponding to K * (892) ± K * 2 (1430) ∓ correlated production.We have observed such correlated production previously in the K + K − π 0 π 0 channel [11]; these results are compared and discussed below.S K 0 L invariant mass distribution for the selected K 0 S K 0 L π + π − events.A clear φ(1020) signal is visible.Fitting with a Gaussian plus polynomial function yields 424 ± 30 φ → K 0 S K 0 L decays, corresponding to about 13% of the events.
We calculate the π + π − invariant mass for events in the φ region, 1.01< m(K 0 S K 0 L ) <1.04 GeV/c 2 , and subtract the non-resonant contribution using events in the sideband 1.04< m(K 0 S K 0 L ) <1.07 GeV/c 2 .We show the resulting m(π + π − ) distribution in Fig. 24(b).It is consistent with those observed in the φπ + π − and φπ 0 π 0 final states [11], where f 0 (980) signals were clearly seen.Fit-ting the m(K 0 S K 0 L ) distribution in bins of the K 0 S K 0 L π + π − mass, we obtain a φπ + π − invariant mass spectrum consistent with those observed in the K + K − π + π − and K + K − π 0 π 0 final states [11].However, the statistical uncertainties are quite large, and so we do not present the distribution or calculate a cross section for this intermediate state.

Final Selection and Backgrounds
This final state contains six charged pions and no neutral particles other than the ISR photon.We consider the events from Sec. III with at least two K 0 S candidates, and the combination of two K 0 S candidates and two charged tracks in each event giving the best χ 2 for a 4C fit under the K 0 S K 0 S π + π − hypothesis (see Sec. V).To reduce the background from multihadronic qq events, we reject events in which both of the charged tracks not in a K 0 S candidate are identified as kaons.
The χ 2 (K 0 S K 0 S π + π − ) distribution for the selected events in the data is shown in Fig. 25 (points), along with that for selected simulated ISR K 0 S K 0 S π + π − events (open histogram), which is normalized to the data in the region χ 2 (K 0 S K 0 S π + π − ) < 10 where the backgrounds and radiative corrections do not exceed 5%.Both distributions are broader than those for a typical 4C χ 2 distribution due to higher-order ISR, and the data include contributions from background processes.

FIG. 25:
The four-constraint χ 2 distributions for K 0 S K 0 S π + π − γ candidate events selected in the data (points) and signal-MC simulation (open histogram) fitted under the K 0 S K 0 S π + π − hypothesis.The cross-hatched histogram represents the simulated background contribution from non-ISR qq events.
The cross-hatched histogram in Fig. 25 represents the background from non-ISR e + e − → qq events.These predominantly contain a hard π 0 , giving a false ISR photon, and have kinematics similar to the signal, giving a peak at low values of χ 2 (K 0 S K 0 S π + π − ).We evaluate this background in a number of E c.m. ranges using the selected data and qq events simulated with JETSET.Combining each ISR photon candidate with all other EMC clusters in the same event, we compare the π 0 signals in the resulting data and simulated γγ invariant mass distributions.The simulation gives an E c.m. dependence consistent with the data, so we normalize its prediction using the overall data-over-MC ratio of π 0 signals, and subtract that from the data.
All remaining background sources are either negligible or yield a χ 2 (K 0 S K 0 S π + π − ) distribution that is nearly uniform over the range shown in Fig. 25.We define signal and control regions, , respectively (see Fig. 25), and use them to estimate and subtract the sum of the remaining backgrounds as described in Sec.VII B. The signal region of Fig. 25 contains 1704 data and 8309 MCsimulated events; the control region contains 219 data and 580 simulated events.
We recalculate the masses of the two K 0 S candidates using the results of the kinematic fit. Figure 26 shows a scatter plot of the invariant mass of one K 0 S candidate versus that of the other for events in the signal region.Any background from events not containing two K 0 S mesons is very low.The m(K 0 S K 0 S π + π − ) distribution for the events in the signal region of Fig. 25 is shown in Fig. 27 as the points.The contributions from non-ISR events and the background estimated from the control region are shown as The K 0 S K 0 S π + π − invariant mass distribution (points) for events in the signal region of Fig. 25.The crosshatched and hatched histograms represent the backgrounds from non-ISR qq events and others estimated from the χ 2 control region of Fig. 25, respectively.The curve represents the smooth empirical fit to the total background used for subtraction.
cross-hatched and hatched histograms, respectively.We fit the sum of all backgrounds with a second order polynomial to reduce fluctuations, and use the result (curve in Fig. 27) for the background subtraction.This gives 1479 signal events with masses between threshold and 4.0 GeV/c 2 .We estimate the systematic uncertainty due to background subtraction to be about 5% of the signal for m(K 0 S K 0 S π + π − ) < 2.5 GeV/c 2 , increasing to about 20% in the 2.5-3.0MeV/c 2 region and 50-70% above 3.0 GeV/c 2 , where background dominates.

B. Simulated detection efficiency
The MC-simulated K 0 S K 0 S π + π − invariant-mass distribution is shown in Fig. 28(a) for events in the signal and control (shaded histogram) regions.The mass dependence of the detection efficiency is shown in Fig. 28(b).The mass interval used, 50 MeV/c 2 per bin, is wider than the 10 MeV/c 2 detector resolution, and the cross section has no sharp structure (except the J/ψ signal, discussed below), so we apply no corrections for the resolution.We apply all the corrections discussed above for data-MC differences in track, K 0 S , and photon detection efficiency.
C. Cross section for e + e − → K 0 S K 0 S π + π − We calculate the e + e − → K 0 S K 0 S π + π − cross section as a function of the effective c.m. energy using Eq. 9 shown in Sec.VII D. The fully corrected cross section is shown in Fig. 29 and listed in Table IV, with statistical uncertainties only.There are no other measurements for this final state.The cross section shows a slow rise from threshold at 1.5 GeV, a maximum value of about 0.5 nb near 2 GeV, and a slow decrease with increasing energy, punctuated by a clear J/ψ signal.The systematic uncertainty is dominated by the uncertainty of the backgrounds, and totals 5% relative at the peak of the cross section, increasing to 20% at 3 GeV, and 50-70% at higher energies.
FIG. 30: The K 0 S π − invariant mass versus the K 0 S π + invariant mass (two entries per event).Figure 30 shows a scatter plot of the K 0 S π − invariant mass versus the K 0 S π + invariant mass, with two entries per event.Clear bands associated with the K * (892) ± are visible here, as are peaks in the projections shown in Fig. 31.The projections also show indications of K * 2 (1430) ± production.
Fitting the projections with a sum of two Breit-Wigner functions and a threshold function yields 827±29 K * (892) + → K 0 S π + and 856±50 K * (892) − → K 0 S π − decays, as well as 116±40 K * 2 (1430) + and 70±34 K * 2 (1430) − decays.The total number of K * (892) ± decays is larger than the number of K 0 S K 0 S π + π − events, indicating correlated production of K * (892) + K * (892) − pairs.We fit the K 0 S π + invariant mass distributon in 0.04 GeV/c 2 bins of the K 0 S π − mass, and show the number of K * (892) + decays in each bin in Fig. 32.A clear K * (892) + signal is observed; a fit yields 742 ± 30 ± 100 The fitted number of K * (892) + events in each 0.04 GeV/c 2 interval of the K 0 S π − mass.The curve represents the result of the fit described in the text, with the hatched area representing the non-resonant component.
These are consistent with the 1:2:1 ratios expected assuming equal production of K 0 S and K 0 L in K * (892) ± decays.
The size of the data sample is not large enough to apply this procedure to every m(K 0 S K 0 S π + π − ) bin and extract the e + e − → K * (892) + K * (892) − cross section.However, considering events with both m(K 0 S π + ) and m(K 0 S π − ) within ±0.15 GeV/c 2 of the nominal K * (892) ± mass [30], we conclude that the K * (892) + K * (892) − contribution almost completely dominates for m(K 0 S K 0 S π + π − ) below 2.5 GeV.For the events outside this box, we show the π + π − and K 0 S K 0 S invariant mass distributions in Fig. 33.The ρ(770) resonance is prominent in the π + π − spectrum, whereas the K 0 S K 0 S spectrum shows no significant structure.The three resonant channels K * (892) + K * (892) − , K * (892) ± K 0 S π ∓ (see Fig. 30), and ρ(770)K 0 S K 0 S dominate the K 0 S K 0 S π + π − cross section within our measured range, and there is a small contribution from K * (1430

Final selection and background
We consider the events from Sec. III with at least two K 0 S candidates, and the combination of two K 0 S candidates and two charged tracks in each event giving the best χ 2 for a 4C fit under the K 0 S K 0 S K + K − hypothesis (see Sec. V).To reduce the background from multi-pionic events, we require that both of the charged tracks not in the K 0 S candidates be identified as kaons.The χ 2 (K 0 S K 0 S K + K − ) distribution for the selected events is shown in Fig. 34 (points) along with that for simulated ISR K 0 S K 0 S K + K − events (open histogram), where the latter distribution is normalized to the data in the region χ 2 (K 0 S K 0 S π + π − ) < 8.There is very little background: simulated ISR events in other channels do not satisfy the selection; there is no significant π 0 peak in the data; and the signal MC describes the data well, even at high χ 2 values.The simulation predicts only a few e + e − → qq → K 0 S K 0 S K + K − π 0 events, which are shown as the hatched histogram in Fig. 34.
We select events with χ 2 (K 0 S K 0 S K + K − ) < 40, obtaining 129 events in the data with masses between threshold and 4.5 GeV/c 2 , and 2544 events in the signal MC simulation.The K 0 S K 0 S K + K − invariant mass distribution is shown as the open histogram in Fig. 35.We do not subtract any background, nor do we assign any systematic uncertainty to account for possible background contributions.The K 0 S K 0 S K + K − invariant mass distribution for data events in the signal region, χ 2 (K 0 S K 0 S K + K − ) < 40 (open histogram).The subset of events with m(K + K − ) < 1.04 GeV/c 2 , predominantly K 0 S K 0 S φ(1020) events, is shown as the shaded histogram.

B. Simulated detection efficiency
The MC-simulated K 0 S K 0 S K + K − invariant-mass distribution is shown in Fig. 36(a) for events in the signal and control (shaded histogram) regions.The mass dependence of the detection efficiency is shown in Fig. 36(b).The mass interval used, 50 MeV/c 2 per bin, is wider than the 10 MeV/c 2 detector resolution, and the cross section has no sharp structure (except the J/ψ signal, discussed below), so we apply no corrections for the resolution.We apply all the corrections discussed above for data-MC differences in track, K 0 S , and photon detection efficiency.FIG.37: The e + e − → K 0 S K 0 S K + K − cross section.Events with invariant mass within 0.05 GeV/c 2 of the J/ψ mass are excluded.

C. Cross section for e
We remove the events within ±0.05 GeV/c 2 of the J/ψ signal (which is discussed below), and calculate the e + e − → K 0 S K 0 S K + K − cross section using Eq. 9.The fully corrected cross section is shown in as a function of energy in Fig. 37 and listed in Table V.There are no previous measurements of this final state.The systematic uncertainties are smaller than the statistical terms and do not exceed 5%.There is also structure for m(K 0 S K 0 S ) near 1.5 GeV/c 2 , which is more visible as a peak in the m(K 0 S K 0 S ) projection of Fig. 38(b).We fit this mass region with a Breit-Wigner plus a second-order polynomial function.An expanded view is shown in Fig. 38(c), along with the result of the fit.We obtain 29±7 events with Breit-Wigner mass and width m = 1.526 ± 0.007 GeV/c 2 Γ = 0.037 ± 0.012 GeV.These parameters may be compared with the averages [30] for the f , and 28.5 ± 5.5 J/ψ → K 0 S K 0 S K + K − decays.Using the respective simulated efficiencies with all the corrections described above, and the differential luminosity, we calculate the products of the J/ψ electronic width and branching fractions to these modes, and list them in Table VI.Using the PDG value of Γ ee (J/ψ ) = 5.55 keV [30], we obtain the corresponding branching fractions, also presented in Table VI.These are the first observations of these J/ψ decay modes and measurements of their branching fractions.They can be compared with B(J/ψ → K + K − π + π − ) = (6.8± 0.3) × 10 −3 [30], which is dominated by the BABAR measurement.
The lines represent the results of the fits described in the text.
A. Internal structure of the J/ψ → K 0 The J/ψ signal in the K 0 S K 0 L π + π − mode has a large non-resonant background (see Fig. 39(a)), and we are unable to quantify the contributions from the K * (892)K 0 S π and φπ + π − intermediate states with reasonable accuracy.The J/ψ → φπ + π − decay rate is relatively well measured [30], dominated by BABAR.
The K 0 S K 0 S π + π − channel has much lower background (see Fig. 39(b)), and we use the 157 events with invariant mass within 30 MeV/c2 of the nominal J/ψ mass to study intermediate states.We use events in the 30 MeV/c 2 intervals on each side of the signal region to estimate a non-J/ψ contribution of 24 events, and to subtract the corresponding contributions from the histrograms that follow.The resulting m(K 0 S π ± ) distribution (four entries per event) is shown in Fig. 40(a).Fitting with two Breit-Wigner (BW) functions plus a polynomial, we obtain 53 ± 14 events containing K * (892)K 0 S π and 35 ± 15 con-taining K * 2 (1430)K 0 S π.To estimate decays to correlated K * (892) + K * (892) − or K * 2 (1430) ∓ K * (892) ± pairs, we consider events from the K * (892) + and K * (892) − bands (see Fig. 30) defined by |m(K 0 S π)−0.892|< 0.15 GeV/c 2 ; a pairing in the overlap region gives only one entry, and there can be as many as two entries per event.Fitting the invariant mass distribution of the other K 0 S π pair, shown in Fig. 40(b), with two BW functions plus a polynomial, we obtain 0.7 ± 5.0 and 8 ± 8 events for the K * (892) + K * (892) − and K * 2 (1430) ∓ K * (892) ± combinations, respectively.Both are consistent with zero, i.e., no correlated production.
For each of these intermediate states we calculate the product of its J/ψ branching fraction, Γ J/ψ ee , and the relevant branching fractions for the intermediate resonances, and list the values in Table VI.Using Γ J/ψ ee = 5.55 eV, known branching fractions [30], and the assumptions that K * mesons decay equally to charged and neutral kaons, and equally to K 0 S and K 0 L (e.g., B K *  S π ± invariant mass distribution (four entries per event) for the K 0 S K 0 S π + π − events under the J/ψ peak with the non-J/ψ contribution subtracted (see text).(b) The distribution of the other m(K 0 S π ∓ ) for those events with one m(K 0 S π ± ) value within 0.15 GeV/c 2 of the K * (892) ± mass (up to two entries per event, pairings in the overlap region taken once).The lines represent the results of the fits described in the text, with the shaded areas representing the combinatorial components.listed in Table VI.The only entry in the PDG tables for any of these channels is B J/ψ→K * (892) + K * (892) − = (1.0 +0.2 −0.4 ) × 10 −3 .Figure 41(a) shows the π + π − invariant mass distribution for the considered events.A clear signal from the ρ(770) resonance is seen, corresponding to J/ψ → ρK 0 S K 0 S decays.The K 0 S K 0 S invariant mass distribution for those events with 0.6 < m(π + π − ) < 1.0 GeV/c 2 , shown in Fig. 41(b), features a narrow spike containing 9.4 ± 4.6 events near 1.53 GeV/c 2 .We observe this same signal when no requirement is placed on the π + π − invariant mass.Attributing this entirely to J/ψ → ρ(770)f ′ 2 (1525) decays, we calculate the measured product and branching fraction, using B(f ′ 2 (1525) → K K) = 0.71 [30], and list them in Table VI.This channel also has no listing in the PDG tables [30].Due to uncertainties in the mass distributions for events without a ρ or f ′ 2 meson, however, we do not attempt to quantify the more inclusive ρK 0 S K 0 S or π + π − f ′ 2 contributions.Using these numbers we calculate the products of Γ J/ψ ee and the relevant branching fractions, and list them in Table VI.Using the PDG values of Γ J/ψ ee , B(φ → K + K − ) = 0.49, and B(f ′ 2 (1525) → K K) = 0.71 [30], we obtain the corresponding branching fractions, also shown in Table VI.Only one value can be compared with an existing PDG listing [30], namly B(J/ψ → f ′ 2 (1525)φ(1020)) = (8 ± 4) × 10 −4 , which has a scale factor of 2.7.Our result can be compared to the MarkII value (4.8 ± 1.8) × 10 −4 , and to the DM2 measurement (12.3 ± 0.26 ± 2.0) × 10 −4 [30].

XII. SUMMARY
We have presented a study of the processes e + e − → K 0 S K 0 L and e + e − → K 0 S K 0 L π + π − at low center-of-mass energies using using events with initial-state radiation (ISR) collected with the BABAR detector.From the dominant e + e − → φγ → K 0 S K 0 L γ process near K 0 S K 0 L threshold, we measure the probability of detecting the K 0 L via its nuclear interaction in the electromagnetic calorimeter with about 0.6% uncertainty, as well as its angular resolution.Using the positions of candidate K 0 L clusters in the calorimeter as input to kinematic fits, we obtain clean samples of K 0 S K 0 L γ and K 0 S K 0 L π + π − γ events, and extract the e + e − → K 0 S K 0 L and e + e − → K 0 S K 0 L π + π − cross sections from threshold to 2.2 and 4 GeV, respectively.
For the K 0 S K 0 L final state, we perform fits to the φ(1020) and φ(1680) resonances, and report the resonance parameters and Γ ee • B(K 0 S K 0 L ) values.The results are consistent with previous measurements, and much more precise for c.m. energies above 1.2 GeV, especially for the φ(1680) mass region.The e + e − → K 0 S K 0 L π + π − cross section is measured for the first time, and is dominated by the K * (892) + K * (892) − intermediate state.Additional contributions from the K * (892) ± K * 2 (1430) ∓ and φπ + π − intermediate states are observed.
We also obtain the first measurements of the e + e − → K 0 S K 0 S π + π − and e + e − → K 0 S K 0 S K + K − cross sections, and provide results from threshold to 4 and 4.5 GeV, respectively.For the former process, we again find the K * (892) + K * (892) − intermediate state to be dominant, and measure a contribution from ρ(770)K 0 S K 0 S .However, no significant contribution from K * (892) ± K * 2 (1430) ∓ is observed.For the latter process, we observe contributions from the K 0 S K 0 S φ(1020) and f ′ 2 (1525)φ(1020) intermediate states.
We observe the J/ψ → K 0 S K 0 L π + π − , K 0 S K 0 S π + π − , and K 0 S K 0 S K + K − decays for the first time, and measure the product of the J/ψ electronic width and branching fraction to each of these modes.

2 FIG. 1 :
FIG.1:The π + π − invariant mass distribution for the selected K 0 S candidates for the data (points) and simulation (histogram).The vertical lines indicate the signal region.

2 FIG. 2 :
FIG. 2:The distribution of constrained recoil mass m c rec , obtained according to Eqs. 1 and 2, for selected γK 0 S candidates.The points represent the data, and the histogram an MC simulation of e + e − → γφ → γK 0 S K 0 L events, normalized to the two most populated bins.

2 FIG. 8 :
FIG.8:The simulated distribution of differences between the reconstructed and generated K 0 S K 0 L invariant mass for the 1 MeV bin of reconstructed mass at the φ peak.

2 FIG. 9 :
FIG.9:The simulated detection efficiency ǫ(s) versus the generated K 0 S K 0 L invariant mass, calculated by dividing the number of events in the signal region of Fig.6by the number generated in each bin.

2 FIG. 10 :
FIG. 10:The K 0 S K 0 L invariant mass distribution in the φ(1020) region.Only statistical uncertainties are shown.The curve represents the result of the fit described in the text.
FIG.11: Two-dimensional plot of the higher cluster energy in a photon-candidate pair versus the corresponding diphoton mass mγγ for all pairs of EMC clusters, containing neither the ISR photon nor the K 0 L candidate.

2 FIG. 13 :
FIG.13:The K 0 S K 0 L invariant mass distribution for data (circles) in the χ 2 control (a) and signal (b) regions (see Fig.12).The histogram in (a) represents MC-simulated signal events in the control region, and that in (b) represents the simulated background from φη, K 0 S K 0 L π 0 , and K 0 S K 0 L 2π 0 events; The squares show the total estimated background, obtained as the difference between the data and simulated distributions in (a), normalized to the data in the signal region.(c) The K 0 S K 0

2 FIG. 14 FIG. 15 :
FIG.14:a)The K 0 S K 0 L mass distributions from the MC simulation for the signal (unshaded) and control (shaded) regions of Fig.12.b) The mass dependence of the simulated net reconstruction and selection efficiency for the K 0 S K 0 L final state in the m(K 0 S K 0 L ) > 1.06 GeV/c 2 region.

2 FIG. 19 :
FIG. 19: (a) The K 0 S K 0 L π + π − invariant mass distribution for MC-simulated signal events in the signal (open histogram) and control region (shaded) of Fig. 18.(b) The net reconstruction efficiency from the simulation.

2 FIG. 23 : 2 FIG. 24 :
FIG.23:The number of K * (892) ± events obtained from fits to the K 0 S π ± invariant mass distribution in each 0.04 GeV/c 2 interval of K 0 L π ∓ mass.The curve represents the result of the fit described in the text, with the hatched areas representing the non-resonant component.

F
Figure24(a)  shows the K 0 S K 0 L invariant mass distribution for the selected K 0 S K 0 L π + π − events.A clear φ(1020) signal is visible.Fitting with a Gaussian plus polynomial function yields 424 ± 30 φ → K 0 S K 0 L decays, corresponding to about 13% of the events.We calculate the π + π − invariant mass for events in the φ region, 1.01< m(K 0 S K 0 L ) <1.04 GeV/c 2 , and subtract the non-resonant contribution using events in the sideband 1.04< m(K 0 S K 0 L ) <1.07 GeV/c 2 .We show the resulting m(π + π − ) distribution in Fig.24(b).It is consistent with those observed in the φπ + π − and φπ 0 π 0 final states[11], where f 0 (980) signals were clearly seen.Fit-

2 FIG. 26 :
FIG.26: Scatter plot of the π + π − invariant mass of one K 0 S candidate versus that of the other K 0 S candidate calculated using the results of the constrained fit.
FIG. 27:The K 0 S K 0 S π + π − invariant mass distribution (points) for events in the signal region of Fig.25.The crosshatched and hatched histograms represent the backgrounds from non-ISR qq events and others estimated from the χ 2 control region of Fig.25, respectively.The curve represents the smooth empirical fit to the total background used for subtraction.
FIG. 28: (a) The K 0 S K 0 S π + π − invariant mass distribution for the MC-simulated signal events in the signal and control (shaded) regions of Fig. 25.(b) The net reconstruction and selection efficiency from the simulation. m.(GeV)

2 FIG. 31 :
FIG.31:The (a) m(K 0 S π + ) and (b) m(K 0 S π − ) projections of Fig.30.The lines and hatched areas represent the resuls of the fits described in the text and their non-K * components, respectively.

χ 2 FIG. 34 :
FIG. 34:The four-constraint χ 2 distributions for K 0 S K 0 S K + K − γ candidate events in the data (points) and signal MC simulation (open histogram) fitted under the K 0 S K 0 S K + K − hypothesis.The cross-hatched histrogram represents the simulated contribution from non-ISR qq events.
FIG. 35:The K 0 S K 0 S K + K − invariant mass distribution for data events in the signal region,χ 2 (K 0 S K 0 S K + K − ) < 40 (open histogram).The subset of events with m(K + K − ) < 1.04 GeV/c 2 , predominantly K 0 S K 0 S φ(1020) events, is shown as the shaded histogram.
FIG. 36: (a) The K 0 S K 0 S K + K − invariant mass distribution for the MC-simulated signal events in the signal and control (shaded) regions of Fig. 34.(b) The net reconstruction and selection efficiency from the simulation.

D
Figure38(a) shows a scatter plot of the K + K − invariant mass versus the K 0 S K 0 S invariant mass for all selected events.A strong φ(1020) band is evident.Requiring m(K + K − ) < 1.04 GeV/c 2 , we obtain the contribution from φK 0 S K 0 S events shown in Fig35as the shaded histogram.This mode dominates at all masses.There is also structure for m(K 0 S K 0 S ) near 1.5 GeV/c 2 , which is more visible as a peak in the m(K 0 S K 0 S ) projection of Fig.38(b).We fit this mass region with a Breit-Wigner plus a second-order polynomial function.An expanded view is shown in Fig.38(c), along with the result of the fit.We obtain 29±7 events with Breit-Wigner mass and width

′ 2 (
Figures39(a), (b), and (c) show expanded views of the mass distributions in Figs.18(c), 27, and 35, respectively, in the J/ψ mass region.Fitting with Gaussian plus second-order polynomial functions yields 248 ± 27J/ψ → K 0 S K 0 L π + π − decays, 133 ± 13 J/ψ → K 0 S K 0 S π + π − decays, and 28.5 ± 5.5 J/ψ → K 0 S K 0 S K + K − decays.Using the respective simulated efficiencies with all the corrections described above, and the differential luminosity, we calculate the products of the J/ψ electronic width and branching fractions to these modes, and list them in TableVI.Using the PDG value of Γ ee (J/ψ ) = 5.55 keV[30], we obtain the corresponding branching fractions, also presented in TableVI.These are the first observations of these J/ψ decay modes and measurements of their branching fractions.They can be compared with B(J/ψ → K + K − π + π − ) = (6.8± 0.3) × 10 −3[30], which is dominated by the BABAR measurement.

2 FIG. 39 :
FIG. 38: (a) TheK + K − versus K 0 S K 0 S invariant mass for all selected K 0 S K 0 S K + K − events in the data.(b)The m(K 0 S K 0 S ) projection of (a).(c) An expanded view of (b) in which the line represents the result of the fit described in the text.
FIG. 40: (a)The K 0 S π ± invariant mass distribution (four entries per event) for the K 0 S K 0 S π + π − events under the J/ψ peak with the non-J/ψ contribution subtracted (see text).(b) The distribution of the other m(K 0 S π ∓ ) for those events with one m(K 0 S π ± ) value within 0.15 GeV/c 2 of the K * (892) ± mass (up to two entries per event, pairings in the overlap region taken once).The lines represent the results of the fits described in the text, with the shaded areas representing the combinatorial components.
FIG. 42: (a) TheK + K − versus K 0 S K 0 S invariant mass for the K 0 S K 0 S K + K − events under the J/ψ peak (see text).(b) The K 0 S K 0 S invariantmass distribution for the events in (a) with m(K + K − ) < 1.04 GeV/c 2 .The line represents the result of the fit described in the text.

Figure 42 ′ 2 (
Figure42(a) shows the K + K − versus K 0 S K 0 S invariant mass for the 30 K 0 S K 0 S K + K − events with total invariant mass within 30 MeV/c 2 of the nominal J/ψ mass, 29±6 of which are J/ψ events.Horizontal and vertical bands are visible, corresponding to the φ(1020) and f ′ 2 (1525) resonances, respectively.We select 20 J/ψ → φ(1020)K 0 S K 0 S candidate decays by requiring m(K + K − ) < 1.04 GeV/c 2 , and plot their m(K 0 S K 0 S ) distribution in Fig.42(b).Fitting with a

TABLE I :
Summary of corrections and systematic uncertainties for the measurement of the e + e − → K 0 S K 0 L process in the φ resonance region.

TABLE II :
Summary of the e + e − → KSKL cross section measurement.Uncertainties are statistical only.

TABLE III :
Summary of the e + e − → KSKLπ + π − cross section measurement.Uncertainties are statistical only.

TABLE IV :
Summary of the e + e − → KSKSπ + π − cross section measurement.Uncertainties are statistical only.

TABLE V :
Summary of the e + e − → KSKSK + K − cross section measurement.Uncertainties are statistical only.

TABLE VI :
Summary of the J/ψ parameters obtained in this analysis.