Observation of the Rare Decay B+ ->K+ pi0 pi0

We report an analysis of charmless hadronic decays of charged B mesons to the final state K+ pi0 pi0, using a data sample of 470.9 +/- 2.8 million BBbar events collected with the BaBar detector at the Upsilon(4S) resonance. We observe an excess of signal events with a significance above 10 standard deviations including systematic uncertainties and measure the branching fraction to be BF(B+ ->K+ pi0 pi0) = (15.5 +/- 1.1 +/- 1.6) x 10^(-6), where the uncertainties are statistical and systematic, respectively.

We report an analysis of charmless hadronic decays of charged B mesons to the final state K + π 0 π 0 , using a data sample of 470.9 ± 2.8 million BB events collected with the BABAR detector at the Υ (4S) resonance. We observe an excess of signal events with a significance above 10 standard deviations including systematic uncertainties and measure the branching fraction to be B B + → K + π 0 π 0 = (15.5 ± 1.1 ± 1.6) × 10 −6 , where the uncertainties are statistical and systematic, respectively.
Recent measurements of rates and asymmetries in B → Kπ decays have have generated considerable interest because of possible hints of new physics contributions [1,2]. Unfortunately, hadronic uncertainties prevent a clear interpretation of these results in terms of physics beyond the Standard Model (SM). A data driven approach, involving measurements of all observables in the B → Kπ system can in principle resolve the theoretical situation, but much more precise measurements (i.e. much larger data samples) will be needed [3][4][5].
It is interesting to study the related decays to pseudoscalar-vector final states B → K * π and B → Kρ [6][7][8]. In Table I we review the existing experimental measurements of the channels in the B → K * π system. It is evident that improved measurements of the K * + π 0 [9] decay are needed.
Due to the non-negligible width of the K * resonances, the quasi-two-body modes are best studied in the analysis of the three-body Dalitz plot. The four K * π decays populate six Kππ Dalitz plots (the four Kρ decays also produce four of the same six final states). Of these, Dalitz plot analyses of K + π + π − [16,17], K 0 S π + π − [14,20] and  [12,19] K + π − π 0 [12] have been performed to date. The first two of these have shown the presence of a poorly-understood structure, dubbed the f X (1300), in the π + π − invariant mass distribution. A study of the K + π 0 π 0 Dalitz plot would help to elucidate the nature of this peak, since even-spin states will populate both Kπ + π − and Kπ 0 π 0 (assuming isospin symmetry), while odd-spin states cannot decay to π 0 π 0 . Knowledge of the K + π 0 π 0 Dalitz plot may also help to clarify the interpretation of the inclusive time-dependent analyses [21] of B 0 → K 0 S π 0 π 0 [22,23]. Currently, these results show the largest deviation, albeit with a large uncertainty, among hadronic b → s penguin-dominated decays [10] from the naïve Standard Model expectation that the time-dependent CP violation parameter should be given by S CP ≈ −η CP sin(2β), where η CP is the CP eigenvalue of the final state (+1 for K 0 S π 0 π 0 ) and β is an angle of the Cabibbo-Kobayashi-Maskawa [24,25] unitarity triangle. Such deviations could be caused by new physics, but in order to rule out the possibility of sizeable corrections to the Standard Model prediction, better understanding of the population of the Kπ 0 π 0 Dalitz plots will be necessary.
The data used in the analysis, collected with the BABAR detector [28] at the PEP-II asymmetric energy e + e − collider at SLAC, consist of an integrated luminosity of 429 fb −1 recorded at the Υ (4S) resonance ("on-peak") and 45 fb −1 collected 40 MeV below the resonance ("offpeak"). The on-peak data sample contains the whole BABAR dataset of 470.9 ± 2.8 million BB events.
We reconstruct B + → K + π 0 π 0 decay candidates by combining a K + candidate with two neutral pion candi-dates. The K + candidates are required to have a minimum transverse momentum of 0.05 GeV/c and to be consistent with having originated from the interaction region. Separation of charged kaons from charged pions is accomplished with energy-loss information from the tracking subdetectors, and the Cherenkov angle and number of photons measured by a ring-imaging Cherenkov detector. The efficiency for kaon selection is approximately 80 % including geometrical acceptance, while the probability of misidentification of pions as kaons is below 5 % up to a laboratory momentum of 4 GeV/c. Neutral pion candidates are formed from pairs of photons with laboratory energies above 0.05 GeV and lateral moments between 0.01 and 0.6. We require that the mass of the reconstructed π 0 is within the range 0.115 GeV/c 2 < m γγ < 0.150 GeV/c 2 and that the absolute value of the cosine of the decay angle in the π 0 rest frame is less than 0.9. We exclude candidates that are consistent with the B + → K 0 S K + , K 0 S → π 0 π 0 decay chain by rejecting events that contain a candidate that satisfies 0.40 GeV/c 2 < m π 0 π 0 < 0.55 GeV/c 2 . This veto has a signal efficiency of at least 96 % for any charmless resonant decay and is almost 100 % efficient for nonresonant B + → K + π 0 π 0 and B + → χ c0 K + decays.
Due to the presence of two neutral pions in the final state, there is a significant probability for signal events to be misreconstructed, due to low momentum particles being exchanged with particles from the decay of the other B meson in the event. We refer to these as "self-crossfeed" (SCF) events, as opposed to correctly reconstructed (CR) events. Using a classification based on Monte Carlo information, we find that in simulated events the SCF fraction depends strongly on the Dalitz plot distribution of the signal, and ranges from 2 % for B + → χ c0 K + decays to 30 % for B + → f 2 (1270)K + decays.
To suppress the dominant background contribution, which arises from continuum e + e − → qq (q = u, d, s, c) events, we employ a neural network that combines four variables commonly used to discriminate jet-like qq events from the more spherical BB events. These are the ratio of the second to the zeroth order momentumweighted angular moment, the absolute value of the cosine of the angle between the B direction and the beam axis, the absolute value of the cosine of the angle between the B thrust axis and the beam axis, and the absolute value of the output of a neural network used for "flavour tagging" (i.e. for distinguishing B from B decays using inclusive properties of the decay of the other B meson in the Υ (4S) → BB decay [29]). The first three quantities are calculated in the center-of-mass (CM) frame. We apply a loose criterion on the neural network output (NN out ) which retains approximately 90 % of the signal while rejecting approximately 82 % of the qq background.
In addition to NN out , we distinguish signal from background events using two kinematic variables: the difference ∆E between the CM energy of the B candi-date and √ s/2, and the beam-energy-substituted mass where √ s is the total CM energy and p B is the momentum of the candidate B meson in the CM frame. The signal m ES distribution for CR events is approximately independent of the B + → K + π 0 π 0 Dalitz plot distribution and peaks near the B mass with a resolution of about 3 MeV/c 2 . We select signal candidates that satisfy 5.260 GeV/c 2 < m ES < 5.286 GeV/c 2 . The CR signal ∆E distribution peaks near zero, but has a resolution that depends on the signal Dalitz plot distribution, which is a priori unknown. To avoid possible biases [30] we apply tighter selection criteria, −0.15 GeV < ∆E < 0.05 GeV, which have an efficiency of about 80 % for signal while retaining only about 30 % of the background (both compared to the looser requirement |∆E| < 0.30 GeV), and do not use ∆E in the fit described below.
The efficiency for signal events to pass all the selection criteria is determined as a function of position in the Dalitz plot. Using a Monte Carlo (MC) simulation in which events uniformly populate the phase-space, we obtain an average efficiency of approximately 16 %, though values as low as 8 % are found near the corners of the Dalitz plot.
An average of 1.3 B candidates is found per selected event. In events with multiple candidates we choose the one with the smallest value of a χ 2 variable formed from the sum of the χ 2 values of the two π 0 candidate masses. This procedure has been found to select the best reconstructed candidate more than 90 % of the time, and does not bias our fit variables.
We study residual background contributions from BB events using MC simulations. It is found that these events can be combined into four categories based on their shapes in m ES and ∆E. The first category comprises two-body modes (mainly B + → K + π 0 ); the second contains three-body modes (mainly B + → K * + γ and B + → π + π 0 π 0 ); the third and fourth are composed of higher multiplicity decays (many possible sources, with or without intermediate charmed states) with missing particles, and are distinguished by the absence or presence of a peak in the m ES distribution respectively. Based on the MC-derived efficiencies, total number of BB events, and known branching fractions [10, 31], we expect 70 ± 9, 39 ± 18, 1090 ± 40 and 170 ± 30 events in the four respective categories.
To obtain the B + → K + π 0 π 0 signal yield, we perform an unbinned extended maximum likelihood fit to the candidate events using two input variables: m ES and NN out . For each component j (signal, qq background, and the four BB background categories), we define a probability density function (PDF) where i denotes the event index. The signal component is further separated into CR and SCF parts where f SCF is the SCF fraction. The extended likelihood function is where n j(k) is the yield of the event category j(k).
For the signal, the m ES PDFs for CR and SCF are described by an asymmetric Gaussian with power-law tails and a third order Chebychev polynomial, respectively. Both CR and SCF NN out PDFs are described by one-dimensional histograms. We fix the shape parameters to the values obtained from the B + → K + π 0 π 0 phase-space MC sample, after adjusting them to account for possible differences between data and MC simulations determined with a high statistics control sample of B + → D 0 ρ + → K + π − π 0 π + π 0 decays. For the continuum background, we use an ARGUS function [32] to parameterize the m ES shape. The continuum NN out shape is modelled with a parametric step function function with 20 bins. One-dimensional histograms are used as nonparametric PDFs to represent all fit variables for the four BB background components. The free parameters of our fit are the yields of signal and continuum background together with the parameters of the continuum m ES and NN out PDFs. All the yields and PDF shapes of the four BB background categories are fixed based on MC simulations.
The results of the fit are highly sensitive to the value of f SCF , which depends strongly on the Dalitz plot distribution of signal events and cannot be determined directly from the fit. To circumvent this problem, we adopt an iterative procedure. We perform a fit with f SCF fixed to an initial value. We then construct the signal Dalitz plot from the signal probabilities for each candidate event ( s Weights) calculated with the s Plot technique [33], and determine the corresponding average value of f SCF . We then fit again with f SCF fixed to the new value, and repeat until the obtained values of the total signal yield (CR + SCF) and f SCF are unchanged between iterations. This method was validated using MC and was found to return values of f SCF that are accurate to within 3 % of the nominal SCF fraction. Convergence is typically obtained within 3 iterations.
We cross-check our analysis procedure using the high statistics control sample described above. We impose selection requirements on the D and ρ candidates' invariant masses: 1.84 GeV/c 2 < m K + π − π 0 < 1.88 GeV/c 2 and 0.65 GeV/c 2 < m π + π 0 < 0.85 GeV/c 2 . We fit the data with a likelihood function that includes components for the control channel, and for backgrounds from BB and qq backgrounds. We find a yield consistent within statistical uncertainties with the expectation based on worldaverage product branching fractions [31].
We apply the fit method described above to the 31 673 selected candidate B + → K + π 0 π 0 events. Convergence is obtained after four iterations with a yield of 1220 ± 85 signal events and a SCF fraction of 9.7 %. The results of the fit are shown in Fig. 1. The statistical significance of the signal yield, given by √ 2∆ ln L where ∆ ln L is the difference between the negative log likelihood obtained assuming zero signal events and that at its minimum, is 15.6 standard deviations (σ). Including systematic uncertainties (discussed below), the significance is above 10 σ.
The B + → K + π 0 π 0 branching fraction is determined from the result of the fit, dividing the signal s Weights by event-by-event efficiencies that take the Dalitz-plot position dependence into account, and summing them to obtain an efficiency-corrected signal yield of 7427 ± 518 events. We further correct for the effect of the K 0 S veto and a bias in the fitted signal yield, as determined from Monte Carlo pseudoexperiments generated with a signal component with the same values of the yield and SCF fraction as found in the fit to data. Finally, we divide by the total number of BB events in the data sample to obtain our measurement of the branching fraction B B + → K + π 0 π 0 = (15.5 ± 1.1 ± 1.6) × 10 −6 , where the first uncertainty is statistical and the second is systematic.
We assign systematic uncertainties due to (i) uncertainties in CR signal m ES PDF shapes (0.8 %) evaluated using the B + → D 0 ρ + → K + π − π 0 π + π 0 control sample; (ii) uncertainties in CR signal and BB background NN out PDF shapes (4.9 %) evaluated using uncertainties in the data/MC ratio determined from the control sample and applying them to all PDFs in a correlated manner; (iii) uncertainties in the SCF signal m ES and NN out PDF shapes (1.7 % and 0.7 %, respectively) evaluated considering a range of SCF shapes corresponding to different signal Dalitz plot distributions; (iv) uncertainty in the SCF fraction (2.5 %) from varying the value used in the fit within a range of uncertainty determined from Monte Carlo pseudoexperiment tests of our iterative fitting procedure; (v) uncertainties in the BB background PDFs due to finite MC statistics (0.8 %), determined by varying the contents of the bins of the histograms used to describe the PDFs within their errors; (vi) uncertainties in the BB background m ES PDF shapes due to data/MC differences (1.6 %), evaluated by smearing the PDFs with a Gaussian with parameters determined from the control sample; (vii) uncertainties in the fixed BB background yields (1.4 %), evaluated by varying these within their uncertainties; (viii) uncertainty in the correction due to fit bias (1.8 %), which corresponds to half the correction combined in quadrature with its error; (ix) uncertainties in the efficiency, with contributions from tracking (0.4 %), kaon identification (1.0 %), neutral pion reconstruction (3.0 % per neutral pion, so 6.0 % in total), ∆E (4.0 %) and NN out (3.0 %) selection requirements, the K 0 S veto correction (2.0 %); (x) uncertainty in the number of BB pairs in the data sample (0.6 %). The total systematic uncertainty on the branching fraction is 10.4 %. Including only systematic uncertainties that affect the yield, the total is 6.5 %. Table II summarizes the systematic  contributions. In summary, using the full BABAR data sample of 429 fb −1 collected at the Υ (4S) resonance, we observe charmless hadronic decays of charged B mesons to the final state K + π 0 π 0 . The signal has a significance above 10 σ, after taking systematic effects into account. We measure the branching fraction to be B B + → K + π 0 π 0 = (15.5 ± 1.1 ± 1.6) × 10 −6 . This is the first step towards understanding the composition of the Dalitz plot of this decay and measuring the properties of contributing quasi-two-body modes.
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