Search for the decay B 0 s → π 0 π 0 at Belle

We report the results of the ﬁrst search for the decay B 0 s → π 0 π 0 using 121 . 4 fb − 1 of data collected at the Υ(5S) resonance with the Belle detector at the KEKB asymmetric-energy e + e − collider. We observe no signal and set a 90% conﬁdence level upper limit of 7 . 7 × 10 − 6 on the B 0 s → π 0 π 0 decay branching fraction.

The study of heavy-flavored hadrons decaying to hadronic final states provides an important input for understanding the interplay between strong and weak interactions.These type of decays involving weak annihilation amplitudes can be a promising place to look for disagreement between theoretical predictions and experimental observations.These decays are highly suppressed and often neglected in theoretical calculations.However, the inclusion of rescattering effects into the theoretical framework naturally enhances their contribution [1].Recently it was observed that the predicted branching fraction for the decay B 0 s → π + π − , which involves topological annihilation diagrams, was substantially smaller than its measured value by the LHCb experiment [2].This dis-crepancy between theoretical prediction and experimental measurement may require some models of strong interaction processes to be revisited [3].In these aspects, searches for decays involving weak annihilation amplitudes become important and necessary.
Within the standard model (SM), the decay B 0 s → π 0 π 0 proceeds via the W -exchange and "penguin" annihilation amplitudes, as shown in Fig. 1.Theoretical calculations based on the Flavor Diagram Approach (FDA) [4], perturbative Quantum Chromodynamics (pQCD) [5] , and QCD factorization [6] predict the branching fraction (B) to be (0.40 ± 0.27) × 10 −6 , (0.28 ± 0.09) × 10 −6 , and (0.13 ± 0.05) × 10 −6 , respectively.The only measurement for this decay was made by the L3 experi-ment in 1995, which reported an upper limit (UL) of B < 2.1 × 10 −4 at 90% confidence level (CL) [7].The search for the decay B 0 s → π 0 π 0 [8] described in this Letter is based on a data sample of 121.4 fb −1 collected at the Υ(5S) resonance using the Belle detector.The Belle detector at the KEKB [9] asymmetric-energy e + e − collider is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector, a 50layer central drift chamber, an array of aerogel threshold Cherenkov counters, a barrel-like arrangement of timeof-flight scintillation counters, and a CsI(Tl) crystalbased electromagnetic calorimeter (ECL) located inside a super-conducting solenoid coil that provides a 1.5 T magnetic field.An iron flux-return outside the coil is instrumented to detect K 0 L mesons and identify muons.A detailed description of the Belle detector can be found elsewhere [10,11].The analysis relies on the ECL component of the detector for the reconstruction of the photons in the B 0 s → π 0 π 0 decay final state.The production cross-section of the Υ(5S) resonance at the e + e − centre of mass (c.m.) energy of 10.86 GeV is σ [12], and the fraction of b b events giving rise to B 0 s production modes, B ( * )0 s B( * )0 s , is measured to be f s = (0.201 ± 0.031) [13].There are three kinematically allowed modes of production of B 0 s mesons: The production fractions from the former two are (87.0 ± 1.7)% and (7.3 ± 1.4)%, respectively [12], while the remaining fraction is from the B 0 s B0 s mode.The B * 0 s decays to B 0 s by radiating a low-energy photon that is usually not identified due to its poor reconstruction efficiency.The number of events with B 0 s B0 s is, therefore, estimated to be • f s = (8.30±1.34)×10 6.We employ a "blind" analysis procedure to leave out the experimenter's biases and develop our analysis strategy with Monte-Carlo (MC) samples.In a "blind" analysis, the signal region is kept hidden until the selection criteria are finalized.The signal MC samples are generated with EvtGen [14] and simulated with GEANT3 [15] to model all possible detector effects.Background studies are performed with MC samples six times larger than the integrated luminosity of data.The analysis procedure is validated with a control sample of B 0 d → π 0 π 0 decays produced at the Υ(4S) resonance, which closely resembles the signal.
We reconstruct B 0 s → π 0 π 0 with π 0 → γγ.Photon candidates are reconstructed from ECL clusters that do not match any charged track and have energy greater than 50 (100) MeV in the ECL's barrel (end-caps) region.The forward end-cap, barrel, and backward endcap regions of the ECL are given by 12 • < θ < 31.4 • , 32.2 • < θ < 128.7 • , and 131.5 • < θ < 157.2 • , respectively, where θ is the polar angle in the laboratory frame with respect to the detector axis, in the direction opposite to the e + beam.To remove the off-time (radiative) Bhabha and e + e − → γγ events, a timing criterion based on the beam collision time is applied, which is determined at the trigger level for each candidate event.The invariant mass of the two-photon combination must lie in the range of 118 MeV/c 2 < m(γγ) < 152 MeV/c 2 , corresponding to ±2.4 standard deviations (σ) of the invariant mass resolution around the nominal π 0 mass [13].A mass-constrained fit is subsequently performed to improve the π 0 momentum resolution.
To further select the B 0 s candidates, we apply selection criteria on their beam-energy-constrained mass , where E beam is the beam energy, p reco and E reco are the momentum and energy, respectively, of the reconstructed B 0 s candidate.The world average value is used for the mass of the B 0 s meson, m B 0 s [13].A B 0 s candidate is retained for further analysis only if it satisfies the requirement that 5.300 GeV/c 2 < M bc < 5.434 GeV/c 2 and −0.60 GeV < ∆E < 0.15 GeV.The backgrounds near the Υ(5S) resonance which can affect the analysis are: continuum (e + e − → q q, q = u, d, s, c), B s → ρ + ρ − and B 0 s → K 0 s π 0 show that their contributions are negligible.We also find no bsbs and nonbsbs background after applying all of the aforementioned selection criteria.Background MC studies, therefore, reveal the dominance of continuum background over the other types of background.Their suppression requires topological variables, which classify the signal and the continuum background based on their event shape variables in the e + e − c.m. frame.
In signal events, B 0 s pairs are produced with small momenta, and the distribution of their decay products tends to be spherical.In contrast, the quark pairs of the continuum background are produced with a significant amount of momentum; therefore, their decay product distribu-tion has a jet-like topology.A neural network algorithm (NN) [16] is employed to suppress the continuum background.The input of the NN includes sixteen modified Fox-Wolfram moments [17], and cos θ T (see Section 9.3 in [18]) to provide additional discrimination between the signal and the continuum background.The angle θ T is defined as the angle between the thrust axis of the signal B 0 s candidate and the thrust axis of the remainder of the events.The NN was trained on MC samples with consistency checks to ensure no over-training.
The choice of the selection criterion on the output of the NN, C NN , is determined based on a Punzi's figureof-merit (FOM) optimization [19], where the significance level is set to three standard deviations.The C NN distributions for the continuum background and the signal lies in the range of [−1, +1], where the continuum backgrounds peak at −1 and the signal candidates at +1.We require C NN to be greater than 0.90 for this analysis.This condition removes 99% of the continuum background with a signal loss of 53%.To facilitate the data modelling, C NN was transformed to another variable, C NN using the following formula where C NN(min) = 0.90 and C NN(max) is the maximum value of C NN obtained from the NN distribution.
After applying the selection criteria described above, 10.3% of signal MC events have more than one candidate.We select the best B 0 s candidate by choosing the one with the smallest sum of the χ 2 of the mass-constrained fits to the two π 0 s.This method chooses the correct B 0 s candidate 56% of the time.The misreconstructed fraction of events after applying all the selection criteria is found to be negligible; hence, they are not treated separately.The overall signal reconstruction efficiency in this analysis is (12.69 ± 0.05)%.
To extract the signal yield, we perform a threedimensional (3D) unbinned extended maximum likelihood (ML) fit to M bc , ∆E , and C NN .The likelihood function is defined as where P j (M bc , ∆E , C NN ) is the PDF of the signal or background component (specified by index j), n j is the yield of this component, i represents the event index, and N is the total number of events in the sample.
The linear correlation coefficients among M bc , ∆E , and C NN are found to be below 3% in the signal region.Consequently, each of the 3D PDFs describing the signal and background contributions are assumed to factorize as P j ≡ P j (M bc )P j (∆E )P j (C NN ).These factorized PDFs are modelled using large signal and background MC sam-ples.The signal M bc PDF consists of three PDFs corresponding to the three B 0 s production channels.Each of them is again separately modelled from large MC samples.They are then combined according to their production fractions [12] to produce the final signal PDF for the M bc variable.The PDF used for parametrizing B 0 s B0 s is a sum of two Gaussian distributions with a common mean, while each of B 0 s B * 0 s or B * 0 B0 s , and B * 0 s B * 0 s , are parametrized using a sum of a Gaussian function and an empirical PDF shape known as the Crystal Ball function [20].The signal ∆E variable, for all the three B 0 s channels, is modelled using the Crystal Ball function, which is modified for this analysis to include the asymmetric nature of the distribution about the mean position.The output from the NN is parametrized using a Gaussian and an asymmetric (bifurcated) Gaussian PDF for the signal C NN variable.Unlike the signal PDF parameters for the M bc variable, which is different for the three B 0 s sources, ∆E and C NN variables take the same parameter values for the three B 0 s production channels.The continuum background distribution of the M bc variable is modelled through an empirically determined parametrized background shape referred to as the ARGUS function [21].The continuum background is parametrized using a first-order Chebychev polynomial and a sum of two Gaussian distributions for the ∆E and C NN variables, respectively.All the signal parameters and the background ARGUS endpoint are fixed to their best fit values obtained from 1D fits to the MC simulated events.In contrast, all other background parameter values and the signal and background yields are floated.The PDFs used for modelling the signal and continuum background are listed in Table I.
TABLE I. PDFs used to model the M bc , ∆E , and C NN distributions.The notations G, BG, CB, ACB, CP, and A correspond to Gaussian, Bifurcated Gaussian, Crystal Ball, Asymmetric Crystal Ball, Chebyshev polynomial, and AR-GUS functions, respectively.

Fit component
To validate our analysis, we use the Belle data sample collected at the Υ(4S) to reconstruct the decay B 0 d → π 0 π 0 by applying similar event selection criteria.The results of the fit to Υ(4S) data are shown in Fig. 2, where each fit projection is plotted after additional selection criteria are applied as described in the caption.We calculate the branching fraction, B(B 0 d → π 0 π 0 ) = (1.18± 0.21) × 10 −6 (where only the statistical uncer-  tainty is shown), which is in good agreement with our previous result [22].
The systematic uncertainties associated with the analysis are summarized in Table II.The systematic uncertainties due to the fit model are determined via ensemble investigations.To carry out an ensemble study, we generate and simulate 500, 000 signal MC events.We randomly select signal events from this sample for different expected signal yields in data.In addition, background MC events are randomly extracted from the background PDFs based on the expected number of background events in the data.This MC sample that now has statistics equivalent to the expected yields in data is amplified by repeating the above procedure a thousand times.We then perform 3D unbinned extended ML fits on these one thousand pseudo-experiments to obtain pull distributions for each of the expected signal yields in data.The average deviation of a constant function fit to the mean of the pull values from the no bias condition is recognized as a fit bias.
We observe a fit bias of −3.3% and assign it as the corresponding systematic uncertainty.The uncertainty due to fixing the parameter values of the PDFs is determined by varying the best fit parameter values within ±1σ of their statistical uncertainties and measuring the deviation of the signal yields in data.We find a fractional systematic uncertainty of +3.5% −5.2% from this source.Apart from fixing the signal PDF parameters and the background PDF's ARGUS endpoint, we have also fixed the fractions of the B 0 s production channels.We vary these fractions within ±1σ of their measured values [12] and repeat the fit.The observed relative variation +5.2% −3.5% of the signal yield is assigned as the systematic uncertainty.The systematic uncertainty of the signal reconstruction efficiency is 0.4% due to the finite number of signal MC events.The systematic uncertainty due to the efficiency of C NN requirement is estimated from the control sample using a parameter, R. It is defined as the ratio between the efficiency of C NN in data and MC.We assign a corresponding systematic uncertainty of ±3% due to the choice of the selection criteria on the NN output.
The systematic uncertainty for the π 0 selection efficiency is determined to be 2.2% per π 0 using the decay τ − → π − π 0 ν τ .Since this uncertainty is completely correlated for the two π 0 s, a total systematic uncertainty of 4.4% is assigned.We assign a fractional systematic uncertainty of 0.03% on the branching fraction of π 0 → γγ [13].The systematic uncertainty due to the b b production cross-section at Υ(5S) resonance, σ b b is estimated to be ±4.7%[12].In addition, the systematic uncertainty due to the three production charmless processes arising from b b events, f s is assumed to be ±15.4% [13].This uncertainty on f s is the dominant systematic uncertainty associated with any B 0 s measurement at Belle.The fit projections obtained from a 3D unbinned extended maximum likelihood fit in the signal regions are shown in Fig. 3.We obtain 5.7 ± 5.8 signal events and 989 ± 32 continuum background events in our fit to the data.The branching fraction is calculated using where N B 0 s B0 s is the number of B 0 s B0 s pairs; rec and N sig yield are the signal selection efficiency obtained from MC simulation and the signal yield obtained from the fit, respectively; and B is the product of the two π 0 -decay branching fractions [13].
Incorporating the signal yield, N sig yield = (5.3), the branching fraction for B 0 s → π 0 π 0 and its product with f s are calculated to be The first uncertainty is statistical, and the second one is systematic.Without significant signal yield, we calculate the UL on the branching fraction using a Bayesian approach.The UL on the branching fraction is estimated by integrating the likelihood function obtained from the maximum likelihood fit procedure from 0% to 90% of the area under the likelihood curve.The systematic uncertainties are incorporated by convolving the likelihood curve with a Gaussian distribution with a mean of zero and width equivalent to the total systematic uncertainty listed in Table II.The UL on the branching fraction, B(B 0 s → π 0 π 0 ) at 90% CL and the product of the branching fraction with f s , f s × B(B 0 s → π 0 π 0 ), is found to be B(B 0 s → π 0 π 0 ) < 7.7 × 10 −6 f s × B(B 0 s → π 0 π 0 ) < 1.5 × 10 −6 The total systematic uncertainties associated with B(B 0 s → π 0 π 0 ) and f s × B(B 0 s → π 0 π 0 ) are +18.1% −18.4% and +9.5% −10.0%, respectively.The results are summarized in Table III.To summarize, we search for the decay B 0 s → π 0 π 0 using the final Belle data sample available at Υ(5S) resonance, which corresponds to an integrated luminosity of 121.4 fb −1 .We do not observe a significant signal yield, and thus set a 90% CL upper limit on the B 0 s → π 0 π 0 branching fraction of 7.7 × 10 −6 .This is the most stringent UL estimated for this decay representing an orderof-magnitude improvement over the previous result [7] by the L3 experiment in 1995.
as bsbs) and B * B * , B * B, B B, B * B * π, B * Bπ, B Bπ and B Bππ (B = B 0 , B + ) decays (referred as non-bsbs).Additional background MC studies on the peaking background of the types B 0

FIG. 2 .
FIG. 2. Signal enhanced projections of M bc (left), ∆E (middle), and C NN (right) for the control sample, B 0 d → π 0 π 0 .Each plot is generated by applying the signal region selection criteria on the two variables other than the plotted variable.The signal regions for the three variables are as follows, 5.2700 GeV/c 2 < M bc < 5.2895 GeV/c 2 , −0.23 GeV < ∆E < 0.15 GeV, and −3.10 < C NN < 7.61.The dark-filled, red (dotted), black (dash-dotted), and blue (solid) color distributions represent the signal, continuum background, rare B 0 d background (backgrounds arising due to b → u transitions) and total fit function, respectively.Points with error bars represent data.

FIG. 3 .
FIG. 3. Signal enhanced projections of M bc (left), ∆E (middle), and C NN (right) for the analysis, B 0 s → π 0 π 0 .Each plot is generated by applying the signal region selection criteria on the two variables other than the plotted variable.The signal regions for the three variables are as follows, 5.395 GeV/c 2 < M bc < 5.434 GeV/c 2 , −0.310 GeV < ∆E < 0.140 GeV, and −3.901 < C NN < 7.451.The dark-filled, red (dotted), and blue (solid) color distributions represent the signal, continuum background and total fit function, respectively.Points with error bars represent data.The peak in the M bc distribution is due to the dominant B 0 s production channel, B * 0 s B * 0 s (87.0%).The other two production channels, B * 0 s B0 s (7.3%) and B 0 s B0 s (5.7%) are present, but suppressed in the plot.

TABLE II .
Summary of systematic uncertainties.