Search for the decay B0 -->K0_S K0_S K0_L

We present the first search for the decay B0 -->K0_S K0_S K0_L using a data sample of 232 million B Bbar pairs. We find no statistically significant evidence for the non-resonant component of this decay. Our central value for the branching fraction, assuming the true Dalitz distribution is uniform and excluding the phi resonance, is B(B0 -->K0_S K0_S K0_L) = (2.4 +2.7 -2.5 +/- 0.6) x 10^{-6} where the errors are statistical and systematic, respectively. We set a single-sided Bayesian upper limit of B(B0 -->K0_S K0_S K0_L)<7.4 x 10^{-6} at 90% confidence level using a uniform prior probability for physical values. Assuming the worst-case true Dalitz distribution, where the signal is entirely in the region of lowest efficiency, the 90% confidence level upper limit is B(B0 -->K0_S K0_S K0_L)<16 x 10^{-6}.

excluding the φ resonance, is B(B 0 → K 0 S K 0 S K 0 L ) = (2.4 +2.7 −2.5 ± 0.6) × 10 −6 where the errors are statistical and systematic, respectively. We set a single-sided Bayesian upper limit of B(B 0 → K 0 S K 0 S K 0 L ) < 7.4 × 10 −6 at 90% confidence level using a uniform prior probability for physical values. Assuming the worst-case true Dalitz distribution, where the signal is entirely in the region of lowest efficiency, the 90% confidence level upper limit is B(B 0 → K 0 S K 0 S K 0 L ) < 16 × 10 −6 .

I. INTRODUCTION
Measurements from the BABAR and Belle experiments have confirmed the Cabibbo-Kobayashi-Maskawa quark mixing matrix (CKM) mechanism [1] as the dominant source of CP violation in flavor-changing weak interactions [2]. In the Standard Model, the mixing-induced, time-dependent CP asymmetry in B 0 decays from b → sqq penguin transitions should be the same as the precisely measured CP asymmetry in charmonium-K 0 decays, namely sin 2β, to within a few percent. New physics from higher mass scales may contribute to the loop in the penguin diagram, which could significantly alter the CP asymmetry in penguin-dominated B decays [3]. Initial CP asymmetry measurements in b → sqq penguin B decays have suggested a possible violation of this test of the Standard Model [4], [5], [6].
The b → sqq penguin decays fall into two categories. If theqq can beūu, a CKM-suppressed tree-level b → u transition can contribute to the decay in addition to the dominant b → sqq penguin. This introduces some uncertainty in the Standard Model prediction of the CP asymmetry, since the b → u and the penguin amplitudes have different weak phases. On the other hand, decays such as B 0 → K 0 S K 0 S K 0 S [5] and B 0 → K 0 S K 0 S K 0 L are purely b → sss penguin transitions and they can only include b → u decay amplitudes through rescattering, thus the Standard Model uncertainty on the predicted CP asymmetry for these decays is generally smaller [7].
It has been noted that three-body B 0 decays of the form B 0 → P P P ′ , where P and P ′ are spin-0 CP eigenstate neutral particles, are CP eigenstates [8], thus the CP asymmetry in B 0 → K 0 S K 0 S K 0 L is not diluted, as is generally the case in three-body B 0 decays. The B 0 → K 0 S K 0 S K 0 L would be a valuable addition to understanding the b → sqq penguin CP asymmetry anomaly, if the branching fraction is large enough. The resonant φK 0 S contribution to this final state, neglecting interference effects, would not yield a sample of signal events large enough to make an interesting CP asymmetry measurement. In the analysis described below, the expected and observed number of K 0 S K 0 S K 0 L events with a K 0 S K 0 L invari- * Also at Laboratoire de Physique Corpusculaire, Clermont-Ferrand, France † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università della Basilicata, Potenza, Italy ant mass in the region of the φ resonance is about eight events. Prior to our search, the non-resonant component of the K 0 S K 0 S K 0 L final state had not been experimentally investigated. This search was motivated by the possibility of discovering a large non-resonant contribution to the K 0 S K 0 S K 0 L final state. Large non-resonant S-wave amplitudes have been found in the Dalitz amplitude anal- . The branching fraction was recently estimated to be [7] .05 −1.96−2.53−0.06 × 10 −6 including the φ resonance. The difference between the upper and lower limits of this estimate is substantial. Another prediction, based on isospin and Bose symmetry [11], gives a somewhat small branching fraction of

II. THE BABAR DETECTOR AND DATASET
The results presented here are based on data collected with the BABAR detector [12] at the PEP-II asymmetric e + e − collider [13] located at the Stanford Linear Accelerator Center. An integrated luminosity of 211 fb −1 , corresponding to 232 × 10 6 BB pairs, was recorded at the Υ (4S) resonance (center-of-mass energy √ s = 10.58 GeV). Charged particles from the e + e − interactions are detected and their momenta are measured using five layers of double-sided silicon microstrip detectors and a 40-layer drift chamber, both operating in the 1.5-T magnetic field of a superconducting solenoid. Photons, electrons, and hadronic showers from K 0 L interactions are identified with a CsI(Tl) electromagnetic calorimeter (EMC). Further charged particle identification is provided by the specific ionization (dE/dx) in the tracking devices and by an internally reflecting ring imaging Cherenkov detector covering the central region. The steel of the magnetic-flux return for the superconducting solenoid is instrumented with resistive plate chambers (instrumented flux return, or IFR), which are used to identify muons and hadronic showers from K 0 L interactions.

A. Daughter candidate selection
We reconstruct K 0 S candidates through the decay K 0 S → π + π − only. We begin by forming all oppositely-charged combinations of reconstructed tracks in the event. The invariant mass is required to be in the range of 490 to 506 MeV/c 2 , assuming the tracks are pions. In addition, the χ 2 probability of the K 0 S vertex fit must be greater than 1% and the transverse component of the decay length is required to be greater than 8 mm. Finally, the angle between the K 0 S flight direction and the K 0 S momentum vector must be less than 0.2 radians. We reconstruct K 0 L candidates by identifying clusters of energy in the EMC and hits in the IFR that are isolated from all charged tracks in the event. The K 0 L candidates based on IFR clusters must have hits in at least two of the 16 to 19 layers in the detector. The K 0 L candidates detected by a hadronic shower in the EMC must have a calorimeter energy of at least 200 MeV, where this energy is from interpreting the calorimeter signal as an electromagnetic shower. Clusters in the EMC that are consistent with photons from π 0 decay are vetoed by two methods. If the K 0 L candidate EMC cluster forms an invariant mass in the range of 100 to 150 MeV/c 2 with another EMC cluster with a calorimeter energy of at least 100 MeV, it is rejected. If the K 0 L candidate cluster has a calorimeter energy greater than 1.0 GeV and two local maxima (or two "bumps") in the spatial distribution of the energy within the cluster, it may be from a high-energy π 0 , where the electromagnetic showers from the two photons are merged into one cluster. If the two bumps within the cluster form an invariant mass greater than 110 MeV/c 2 , the candidate is rejected.
Further background rejection for EMC K 0 L candidates is achieved by using a neural network trained on signal, qq continuum (where q = u, d, s, c), and BB Monte Carlo samples to distinguish K 0 L clusters from fake clusters. The neural network inputs are the EMC cluster energy and the following six shower-shape variables: • The lateral moment LAT of the shower energy deposition [14] defined as LAT = where the n crystals in the EMC cluster are ranked in order of deposited energy (E i ) , r 0 = 5 cm is the average distance between crystal centers, and r i is the radial distance of crystal i from the cluster center.
• The second radial moment of the shower energy deposition, defined as i E i r 2 i where r i is the radial distance of crystal i from the cluster center.
• The energy sum of a 3×3 block of crystals, centered on the single crystal with the most energy, divided by the larger 5 × 5 block, also centered in the same way.
• The energy of the most energetic crystal in the cluster divided by the energy sum of the 3 × 3 crystal block with the most energetic crystal in the center.
The Zernike moment A n,m is defined as where r i and φ i are the radial and angular separation of crystal i with respect to the cluster center, E tot is the total cluster energy, and R 0 is a cutoff radius of 15 cm. Only K 0 L candidates that pass an optimized cut on the neural network output are retained. This cut has a signal efficiency of 85% and rejects 70% of the EMC fake K 0 L background.

B. B 0 candidate selection
We reconstruct B 0 candidates from selected K 0 L clusters and pairs of K 0 S candidates that do not share any tracks. We require the sum of the K 0 S momentum magnitudes in the center-of-mass frame to be at least 2.1 GeV/c, which ensures consistency with the kinematics of a B 0 → KKK decay. Only the direction of the K 0 L is reconstructed, from the vector defined by the primary vertex and the center of the neutral cluster. We compute the K 0 L momentum by constraining the K 0 We reject fake K 0 L candidates by using the difference between the K 0 L transverse momentum, computed from the B 0 mass constraint, and the transverse momentum along the K 0 L direction that is missing from the event. The reconstructed event missing transverse momentum is calculated without using the K 0 L cluster in the momentum sum and projected along the computed direction of the K 0 L . This missing transverse momentum difference (reconstructed minus calculated) is required to be greater than −0.5 GeV/c. This requirement and the previously mentioned requirement on the EMC K 0 L neural network output were simultaneously optimized to give the greatest signal significance S/ √ S + B, where S and B are the expected number of signal and background events, assuming a signal branching fraction of 5 × 10 −6 [7].
The difference in energy ∆E between the reconstructed B 0 candidate and the beam energy in the center-of-mass frame is the main variable used to distinguish properly reconstructed signal events from combinatoric background. The missing transverse momentum difference distributions for the signal Monte Carlo sample and the background-dominated ∆E sideband (∆E > 0.010 GeV in data) are shown in Figure 1.
We distinguish the non-resonant three-body B 0 decay from two-body B 0 decays to the same final state in this analysis. There are four two-body B 0 decays through charmonium that can go to the final state K 0 are particularly unwanted, since they are from a colorsuppressed tree decay amplitude, not the b → sqq penguin decay amplitude. The χ c modes have unknown B 0 branching fractions. We veto B 0 candidates consistent mass in the range 3.400 to 3.429 GeV/c 2 or 3.540 to 3.585 GeV/c 2 respectively. The combined B 0 and daughter branching fractions for the J/ψK 0 S and ψ(2S)K 0 S modes are 0.062 and 0.016 times 10 −6 respectively. We expect about two events from two-body B 0 decays through charmonium in our final sample.
We also remove B 0 candidates consistent with B 0 → φK 0 S with φ → K 0 S K 0 L by requiring the invariant mass of both K 0 S K 0 L combinations to be above 1.049 GeV/c 2 , though as a cross check, we measure the branching fraction in only the φK 0 S region, where one K 0 S K 0 L combination has an invariant mass less than 1.049 GeV/c 2 . These two-body B 0 decay vetos are 96.6% efficient for the B 0 → K 0 S K 0 S K 0 L signal Monte Carlo sample generated with a uniform true Dalitz distribution.

C. Event selection
The main source of background is from random K 0 S K 0 S K 0 L combinations from continuum qq events. We combine 3 variables in a neural network that is trained using Monte Carlo samples to distinguish signal from continuum events. The first input is the cosine of the angle of the B 0 momentum with respect to the beam axis in the center-of-mass frame θ * B , which is flat for continuum background whereas the signal probability is proportional to sin 2 θ * B . The other two inputs are topological variables that are commonly used to distinguish jet-like continuum events from the more isotropic particle distributions in BB events. The first is the cosine where θ i is the angle with respect to the B 0 thrust axis of track or neutral cluster i, p i is its momentum, and the sum excludes the B 0 candidate daughters. Distributions of the neural network output N N for the signal Monte Carlo sample and the background-dominated ∆E sideband are shown in Figure 2. We require N N > 0.5 to remove events that have little probability of being signal.
After all of the requirements stated thus far, if an event has more than one B 0 candidate, we choose the best one by selecting on the quality of the K 0 L cluster. If there are one or more EMC K 0 L candidates, the best K 0 L candidate is the one with the highest cluster calorimeter energy. If there are no EMC K 0 L candidates, the best K 0 L candidate is the one with the highest number of IFR layers with hits in the K 0 L cluster. If there is more than one K 0 S pair that uses the same (best) K 0 L cluster, the best B 0 candidate is the one with the lowest K 0 S mass χ 2 , defined below where δm i is the difference between the reconstructed invariant mass of K 0 S candidate i and the known K 0 S mass [16] and σ is the invariant mass resolution (2.9 MeV/c 2 ).
Our We use a sample of B 0 → J/ψK 0 L decays from the data to calibrate the reconstruction and selection efficiency and the ∆E resolution. The J/ψ is reconstructed in the e + e − and µ + µ − channels. We apply the same K 0 L selection, projected missing transverse momentum difference, and N N criteria as for our K 0 S K 0 S K 0 L selection. The ratio of the K 0 S K 0 S K 0 L and J/ψK 0 L Monte Carlo efficiencies for these criteria is ǫ K 0 We compare the number of fitted J/ψK 0 L events to the predicted number of events, based on the known branching fractions and the Monte Carlo efficiency. The J/ψK 0 L reconstruction efficiency for all selection criteria is 12.5%. We find N obs = 1420 ± 56 J/ψK 0 L events, consistent with the predicted yield of N 0 = 1426 ± 86. We correct the K 0 S K 0 S K 0 L efficiency by multiplying the MC efficiency by a correction factor of 0.96 ± 0.08, which is defined as The uncertainty on the correction factor includes the uncertainties on the relevant branching fractions for J/ψK 0 L , the statistical error from the J/ψK 0 L fit, and our estimated uncertainty on relating the J/ψK 0 L and K 0 S K 0 S K 0 L selection and reconstruction efficiencies.
The ∆E resolution is better for EMC candidates because the position of the hadronic shower is more precisely measured. We also use the J/ψK 0 L to measure the signal ∆E resolution separately for K 0 L candidates reconstructed in the EMC and IFR. Figure 5 shows the fitted ∆E distributions for EMC and IFR K 0 L candidates.

E. Signal yield determination
We use an extended unbinned maximum likelihood fit to determine the number of signal events in our final sample. The likelihood for an event is the product of probability density functions (PDFs) for the two main discriminating variables in the analysis: ∆E and N N . Separate PDFs are used for EMC and IFR candidates due to the difference in ∆E resolution and the fact that some background channels mostly produce fake signal for EMC K 0 L candidates only. The EMC and IFR K 0 L samples are fitted simultaneously and the relative fraction of EMC signal events is taken from the J/ψK 0 L control sample. Separate PDFs are also used for signal, combinatoric background, and three classes of background from B decays similar to our signal ("peaking backgrounds" described below in Section III F). The combinatoric background includes random K 0 S K 0 S K 0 L combinations from BB events as well as from continuum events. Figure 6 shows the ∆E and N N PDFs for all five fit components. The numbers of signal and combinatoric background events are determined from the fit. The peaking background yields are fixed to values based on known and estimated branching fractions and then varied in the evaluation of systematic uncertainties.
The functional form of the signal PDF for ∆E is a triple Gaussian distribution. The mean and width of the core Gaussian distribution are determined separately for EMC and IFR K 0 L candidates from the J/ψK 0 L control sample and held fixed in the K 0 S K 0 S K 0 L fit. The remaining Gaussian parameters for the signal ∆E PDF are determined from the signal Monte Carlo sample and held fixed. The signal PDF for N N is a 4th-order polynomial. The polynomial coefficients are determined from the signal Monte Carlo sample and held fixed in the fit. The combinatoric background ∆E PDF is an ARGUS function [17]. The N N PDF for combinatoric background is the sum of a 1st-order polynomial and an ARGUS function. The PDF shape parameters for the combinatoric background component are free in the fit.

F. Backgrounds from other B decays
Backgrounds from B decays to final states similar to K 0 S K 0 S K 0 L can look similar to our signal in our discriminating variables ∆E and N N . We call events from these decays "peaking backgrounds." The largest single source of peaking background is the decay B 0 → K 0 S K 0 S K 0 S . One of the K 0 S can decay to π 0 π 0 , where one or more of the photons from either π 0 fakes the K 0 L cluster in the EMC. The world average K 0 S K 0 S K 0 S branching fraction from the Heavy Flavor Averaging Group (HFAG) [6], [5] is B(B 0 → K 0 S K 0 S K 0 S ) = (6.2 ± 0.9) × 10 −6 . We expect 20 K 0 S K 0 S K 0 S events in our fit sample and vary this number by ±5 events in the evaluation of the systematic errors.
Decays of the type KKK * can produce a K 0 S K 0 S K 0 L π final state, which can look similar to our signal if the momentum of the additional π is low. The branching fractions for the relevant KKK * have not been measured. We assume they are each of the same order as our expected signal, based on comparing the relative branching ratios of B → Kπ decays with B → K * π and B → Kρ decays. We estimate a total branching fraction of 30 × 10 −6 for this inclusive final state and vary this assumption by ±15 × 10 −6 in the evaluation of systematic errors. This gives us an expected 70 ± 35 events from KKK * in the fit sample. In the evaluation of the 90% C. L. upper limit on the branching fraction, we set the KKK * event yield to zero. This is the most conservative assumption, since the signal and KKK * event yields are anti-correlated in the fit.
Finally, the decay B 0 → K 0 S K 0 S π 0 can also mimic our signal if the π 0 is misidentified as a K 0 L in the EMC. The branching fraction for this mode has also not been measured. We estimate the branching fraction, relative to our signal, by assuming that the tree-to-penguin ratio for the three-body decays is the same as for B 0 → π 0 π 0 vs B 0 → K 0 π 0 and that our signal branching fraction is about 5 × 10 −6 [7]. Our estimate for the branching fraction is B(B 0 → K 0 S K 0 S π 0 ) ≈ 1.6 × 10 −6 and we vary this assumed branching fraction in the range [0.1, 5.0] × 10 −6 in the evaluation of systematic errors. This gives us an expected 2 +4 −2 events from K 0 S K 0 S π 0 in the fit sample. The ∆E and N N PDFs for the three peaking background components above are determined from Monte Carlo samples. In general, the ∆E shape peaks near ∆E = 0, though there is a substantial tail that extends to high ∆E values. We use the sum of an ARGUS function and a Landau function for the peaking background ∆E PDF. The N N shape is quite similar to the signal shape for all three peaking backgrounds. The form of the N N PDF is a polynomial. The ∆E and N N PDFs for the three peaking background components are shown in Figure 6: c, d, and e. Table I lists the results of the maximum likelihood fit. We find 23 +23 −22 ± 6 events for the K 0 S K 0 S K 0 L signal yield, where the first errors are statistical and the second error is systematic. The maximum likelihood fit bias of +0.3 signal events was evaluated from an ensemble of data sets composed of fully simulated signal and B background Monte Carlo events. The combinatoric background for these datasets was generated from the PDF parameters of the fit to the data. This average bias of 0.3 events and the expected 2.1 events from B 0 decays through charmonium are subtracted from the fitted signal yield.

IV. RESULTS
The systematic errors on the fitted signal yield and the signal branching fraction are listed in Table II  fixed PDF parameters (2.6 events), our uncertainty on the fit bias correction (0.2 events), and the uncertainty on the charmonium background subtraction (1.1 events). The multiplicative contributions are from the uncertainties on the K 0 S and K 0 L reconstruction efficiencies (6.8 and 8.0%), the number of BB events in the sample (1.1%), and the K 0 S → π + π − branching fraction. Figure 7 shows the fitted distributions of ∆E and N N , where a cut has been made on the fit variable not shown (e.g. there is a N N cut applied for the ∆E plot) to enhance the signal. The signal efficiency of this cut is 57% for the ∆E distribution and 71% for the N N distribution. The true distribution of K 0 S K 0 S K 0 L events in the Dalitz plot is unknown. We use the average signal efficiency assuming a uniform true Dalitz distribution (8.1 %) to compute the K 0 S K 0 S K 0 L branching fraction. We found no significant dependence of the signal ∆E resolution or the signal N N shape on the Dalitz plot variables. With these assumptions, we find a K 0 where the first error is statistical and the second error is systematic. The dominant systematic error is from the uncertainties of peaking B backgrounds (23% relative). Figure 8 shows a scan of the negative log likelihood as a function of the K 0 S K 0 S K 0 L branching fraction, where the minimum negative log likelihood has been subtracted. In order to remove any dependence on our estimate of the total KKK * branching fraction, we conservatively fix the KKK * yield to zero in the scan of the log likelihood. We compute a one-sided Bayesian 90% confidencelevel upper limit on the branching fraction assuming a uniform prior probability for positive (physical) branching fraction values. Systematic errors are included by convolving the fit likelihood with a Gaussian distribution with a width corresponding to the total systematic error, excluding the uncertainty on the KKK * yield, since it is fixed to zero in the likelihood scan. The result for the non-resonant three-body branching fraction is Assuming the worst-case true Dalitz distribution, where the signal is entirely in the region of lowest efficiency (4%), the 90% C. L. upper limit on the branching fraction is 16 × 10 −6 .
FIG. 8: Scan of the fit negative log likelihood as a function of the K 0 S K 0 S K 0 L branching fraction, where the minimum negative log likelihood has been subtracted. A uniform true Dalitz distribution has been assumed for the K 0 S K 0 S K 0 L signal. The solid curve includes the systematic uncertainty with yield of KKK * fixed at 0 to calcuate upper limit.
As a cross check on the analysis, we have performed the fit only in the φK 0 S region of the Dalitz plot where the invariant mass of one of the K 0 S K 0 L pairs is less than 1.049 GeV/c 2 . The results are given in the second column of Table I. We find 8.3 +5.5 −4.5 signal events, which corresponds to a branching fraction of B(B 0 → φK 0 S ) = (4.0 +2.6 −2.2 ) × 10 −6 , where the errors are statistical only. This is consistent with the world average value of 0.5 · B(B 0 → φK 0 ) = (4.3 +0.7 −0.6 ) × 10 −6 [16]. We checked our estimation of the KKK * peaking background yield by allowing it to float in the fit. This fit gave a KKK * yield of −54 ± 170 events, which is consistent with our estimation of 70 ± 35 events.

V. SUMMARY
We have searched for the decay B 0 → K 0 S K 0 S K 0 L using 232 million BB pairs recorded by the BABAR experiment. We find no significant evidence for this decay. The central value for the branching fraction, assuming a uniform true Dalitz distribution for the signal and excluding the φ resonance, is B(B 0 → K 0 S K 0 S K 0 L ) = (2.4 +2.7 −2.5 ± 0.6) × 10 −6 , where the first error is statistical and the second I: Results of the fit for yields, branching fraction calculation, and branching fraction upper limit (UL) calculation. The K 0 S K 0 S K 0 L efficiency below assumes a uniform true Dalitz distribution. The maximum likelihood fit bias and charmonium background are subtracted from the signal yield branching fraction calculation.
Assuming the worst-case true Dalitz distribution, where the signal is entirely in the region of lowest efficiency, the upper limit on the branching fraction is 16 × 10 −6 . Our results show that the K 0 S K 0 S K 0 L channel will be of limited use in understanding the b → sqq penguin CP anomaly, due to the low efficiency times branching fraction, which limits the yield of signal events.

VI. ACKNOWLEDGMENTS
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy II: Estimates of systematic errors. Multiplicative systematic errors are in percent while additive systematic errors are in events. The fit yield systematic error is due to fixed fitting parameters, mainly from K 0 S K 0 S K 0 L ∆E PDF parameters, and the uncertainties of fixed peaking B background yields. The fit bias error is one half of the bias. The ccK S/L error is the uncertainty of charmonium background subtraction, estimated as one half of the subtraction.