Measurement of Branching Fraction and Search for $CP$ Violation in $B\to \phi \phi K$

We report the measurement of branching fractions and $CP$-violation asymmetries in $B\to \phi \phi K$ decays based on a $711\,{\rm fb}^{-1}$ data sample containing $772\times 10^6$ $B\bar{B}$ events. The data were recorded at the $\Upsilon (4S)$ resonance with the Belle detector at the KEKB asymmetric-energy $e^+ e^-$ collider. For $B^+ \to \phi \phi K^+$, the branching fraction and $CP$-violation asymmetry measured below the $\eta_{c}$ threshold ($m_{\phi\phi}<2.85\,{\rm GeV}/c^2$) are $[3.43^{\,+\,0.48}_{\,-\,0.46}({\rm stat})\pm 0.22({\rm syst})] \times10^{-6}$ and $-0.02\pm0.11({\rm stat})\pm0.01({\rm syst})$, respectively. Similarly, the branching fraction obtained for $B^0 \to\phi\phi K^0$ below the $\eta_{c}$ threshold is $[3.02^{\,+\,0.75}_{\,-\,0.66} ({\rm stat})\pm \,0.20({\rm syst})]\times10^{-6}$. We also measure the $CP$-violation asymmetry for $B^+ \to\phi\phi K^+$ within the $\eta_{c}$ region ($m_{\phi\phi}\in [2.94,3.02]\,{\rm GeV}/c^2$) to be $+0.12\pm0.12({\rm stat})\pm0.01({\rm syst})$.

PACS numbers: 13.25.Hw, 14.40.Nd B-meson decays to three-body φφK final states proceed via a b → sss loop (penguin) transition, which requires the creation of an additional ss pair. The same final state can also originate from the tree-level process B → η c (→ φφ)K. Figure 1 shows the dominant Feynman diagrams that contribute to these decays. The interference between penguin and tree amplitudes is maximal when the φφ invariant mass lies close to the η c mass (m φφ ∈ [2.94, 3.02] GeV/c 2 ). No CP violation is expected from this interference, as the relative weak phase between the two amplitudes is arg(V tb V * ts /V cb V * cs ) ≈ 0, where V ij denote CKM matrix elements [1]. A potential new physics (NP) contribution to the loop, however, can introduce a nonzero CP -violating phase. In particular, the CP asymmetry can be as large as 40% in the presence of NP [2]. Thus, an observation of large CP violation in B → φφK would indicate the presence of physics beyond the Standard Model. In addition to being an NP probe, the decay is sensitive to the possible production of a glueball candidate near 2.3 GeV/c 2 that can subsequently decay to φφ [3]. We can also search for a structure at 2.35 GeV/c 2 observed in the m φφ distribution in two-photon collisions [4] and dubbed the X(2350). Based on a 78 fb −1 data sample, Belle reported the first evidence for the decay with a branching fraction B(B + → φφK + ) = [2.6 + 1.1 − 0.9 (stat) ± 0.3(syst)] × 10 −6 [5] below the η c threshold (m φφ < 2.85 GeV/c 2 ) [6]. The result was consistent with the corresponding theory prediction, which lies in the range (1.3-4.2)×10 −6 [7,8]. The BaBar experiment performed a measurement of this decay using their full dataset of 464×10 6 BB events [9]. The branching fraction obtained with the same m φφ requirement was B(B + → φφK + ) = (5.6 ± 0.5 ± 0.3) × 10 −6 , about three standard deviations above Belle's result and larger than theoretical estimates. The B 0 → φφK 0 channel was ob-served with a branching fraction of (4.5±0.8±0.3)×10 −6 . BaBar also reported CP asymmetries for charged B decays as −0.10 ± 0.08 ± 0.02 below the η c threshold and +0.09 ± 0.10 ± 0.02 within the η c region.
In this paper, we update our earlier result [5] with a significantly larger data sample containing 772 × 10 6 BB events. The data were collected at the Υ (4S) resonance with the Belle detector [10] at the KEKB asymmetricenergy e + e − collider [11]. The subdetectors relevant for our study are a silicon vertex detector (SVD), a central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), and time-of-flight scintillation counters (TOF). All these are located inside a 1.5 T axial magnetic field.
To reconstruct B + → φφK + and B 0 → φφK 0 decay candidates, we combine a pair of φ mesons with a charged kaon and K 0 S , respectively. All charged tracks except for those from the K 0 S must have a distance of closest approach with respect to the interaction point (IP) of less than 0.2 cm in the transverse r-φ plane, and less than 5.0 cm along the z axis. The z axis is defined as the direction opposite that of the e + beam. We identify charged kaons based on a likelihood ratio R K/π = L K /(L K +L π ), where L K and L π denote the individual likelihood for kaons and pions, respectively. These are calculated using specific ionization in the CDC and information from the ACC and the TOF. A requirement R K/π > 0.6 is applied to select kaon candidates. The kaon identification efficiency, averaged over the momentum range, is 90%, with a pion misidentification rate of about 10%.
We reconstruct the φ candidates from pairs of oppositely charged kaons with an invariant mass in the range 1.00-1.04 GeV/c 2 , corresponding to ±5σ (σ is the width of the mass distribution) around the nominal φ mass [12]. This is referred to as the M KK signal region in the following discussion. The K 0 S candidates are reconstructed from two oppositely charged tracks, assumed to be pions, and are further required to satisfy a criterion on the output of a neural network (NN) algorithm [13]. The algorithm uses the following input variables: the K 0 S momentum in the lab frame; the distance of closest approach along the z axis between the two pion tracks; the flight length in the r-φ plane; the angle between the K 0 S momentum and the vector joining the IP to the K 0 S decay vertex; the angle between the K 0 S momentum in the lab frame and the pion momentum in the K 0 S rest frame; the distances of closest approach in the r-φ plane between the IP and the two pion tracks; the number of CDC hits for each pion track; and the presence or absence of SVD hits for each pion track. We require that the invariant mass lie between 491 MeV/c 2 and 504 MeV/c 2 , which corresponds to a ±3σ window in resolution around the nominal K 0 S mass [12]. B-meson candidates are identified with two kinematic variables: the beam-energy-constrained mass M bc ≡ is the beam energy, and p i and E i are the momentum and energy, respectively, of the i-th decay product of the B candidate. All these quantities are evaluated in the e + e − center-of-mass (CM) frame. We perform a fit for each B candidate, constraining its decay products to originate from a common vertex. Candidate events with M bc ∈ [5.230, 5.289] GeV/c 2 and |∆E| < 0.1 GeV are retained for further study. The M bc requirement corresponds to approximately (−16σ, +3σ) in resolution around the nominal B mass [12], and the ∆E requirement denotes a ±10σ window around zero. We apply such loose requirements on M bc and ∆E as these are used in a maximum-likelihood fit to obtain the signal yield (described later). We define a signal region as M bc ∈ [5.272, 5.289] GeV/c 2 and |∆E| < 0.05 GeV.
After application of the above selection criteria, the average number of B candidates found per event selected in data are 1.7 and 1.6 for B + → φφK + and B 0 → φφK 0 , respectively. In the case of multiple B candidates, we choose the candidate with the lowest χ 2 value for the aforementioned B-vertex fit. From Monte Carlo (MC) simulation the best candidate selection method is found to have an efficiency of 68% (65%) to correctly identify the B-meson candidate in B + → φφK + (B 0 → φφK 0 ) decays. In only about 6% of the total signal events, the B candidate is misreconstructed due to swapping of kaons between the two φ candidates, or of one daughter track with that from the rest of the event. Such misreconstructed events are treated as a part of the signal.
The dominant background is from the e + e − → qq (q = u, d, s, c) continuum process. To suppress this background, observables based on event topology are used. The event shape in the CM frame is expected to be spherical for BB events and jet-like for continuum events. We use an NN [13] to combine the following six variables: a Fisher discriminant formed out of 16 modified Fox-Wolfram moments [14]; the cosine of the angle between the B momentum and the z axis; the cosine of the angle between the B thrust axis [15] and the z axis; the cosine of the angle between the thrust axis of the B candidate and that of the rest of the event; the ratio of the second-to the zeroth-order Fox-Wolfram moments (all quantities are calculated in the CM frame); and the vertex separation along the z axis between the B candidate and the remaining tracks. The NN training and validation are performed with signal and qq MC simulated events. The signal sample is generated with the EvtGen program [16], assuming a uniform distribution over the three-body phase space of the final state.
The neural network output (O NN ) ranges between −1.0 and 1.0, where events near −1.0 (1.0) are more continuum-(signal-) like. We apply a loose criterion O NN > −0.5 to reduce the continuum background. The relative signal efficiency loss due to this requirement is about 6% (3%) for B + → φφK + (B 0 → φφK 0 ) decays, whereas the fraction of continuum events rejected is 76% (66%). As the remainder of the O NN distribution strongly peaks near 1.0 for signal, it is difficult to model with an analytic function. However, the transformed variable where O NN,min = −0.5 and O NN,max ≃ 1.0, has a Gaussian-like distribution that is easier to model. Thus, we use this transformed variable in our signal fit. Backgrounds due to B decays, mediated by the dominant b → c transition, are studied with MC samples of such decays. For both B + → φφK + and B 0 → φφK 0 channels, the M bc and ∆E distributions are found to peak in the signal region. To investigate the source of these contributions, we inspect the m φφ distribution, which displays several peaks corresponding to the η c and other charmonium resonances. To suppress these peaking backgrounds, we exclude candidates for which the m φφ value is greater than 2.85 GeV/c 2 . This requirement also allows us to compare our results with the earlier ones from Belle [5] and BaBar [9]. We calculate the detection efficiencies for candidate events below the η c threshold to be 12.4% and 12.0% for B + → φφK + and B 0 → φφK 0 , respectively.
Charmless backgrounds that do not produce only kaons in the final state may still contribute to the M bc -∆E signal region when a final-state particle is misidentified. These are studied with a BB MC sample in which one of the B mesons decays via b → u, d, s transitions with known or estimated branching fractions [12]. Only 40 events survive from an MC sample equivalent to 50 times the size of the data sample. This small component is combined with the events surviving from b → c transitions to form an overall BB background component. In addition to this BB background that does not peak in M bc or ∆E, we can have contributions from B → φKKK and B → KKKKK decays (described later), which have the same final-state particles as the signal.
The signal yield is obtained with an unbinned extended maximum-likelihood fit to the three variables M bc , ∆E, and O ′ NN . We define a probability density function (PDF) for each event category, i.e., signal, qq, and BB backgrounds: where i denotes the event index, q i is the charge of the B candidate (q i = ±1 for B ± ), and P j and A CP,j are the PDF and CP asymmetry, respectively, for the event category j. The latter is defined as where N B + (N B − ) is the number of B + (B − ) events. We find equal detection efficiencies for the B + (12.3 ± 0.1%) and B − (12.4 ± 0.1%) decays. For neutral B decays, we replace the factor 1 2 (1 − q i A CP,j ) by 1 in Eq. (2). We also do not perform a CP -violation study in this case, since we would need to tag the recoiling B candidate for that, causing further loss in efficiency on top of the small signal yield. As the correlations among M bc , ∆E, and O ′ NN are found to be small ( 5%), the product of three individual PDFs is a good approximation for the total PDF. The extended likelihood function is where n j is the yield of event category j, and N is the total number of candidate events. From the fitted signal yield (n sig ), we calculate the branching fraction as where ε and N BB are the detection efficiency and the number of BB events, respectively. In case of B 0 → φφK 0 , we multiply the denominator by a factor of 1 2 to account for K 0 → K 0 S , as well as by the subdecay branching fraction B(K 0 S → π + π − ) [12]. As the expected yield of the nonpeaking BB background is small, and it is distributed similarly to qq in M bc and ∆E, we merge qq and BB backgrounds into a single component. We find that the difference in the O ′ NN distribution between the two backgrounds contributes a negligible systematic uncertainty. Table I lists the PDF shapes used to model M bc , ∆E, and O ′ NN distributions for various event categories of B → φφK candidates. The yield and PDF shape parameters of the combined background are floated in B + → φφK + . For the neutral channel, however, the background PDF shapes are fixed to their MC values after correcting for small differences between data and simulation, as obtained from the charged decay. Similarly, for the signal components, we fix the M bc , ∆E, and O ′ NN shapes to MC values and correct for small data-MC differences according to values obtained from a control sample of B + → D + s D 0 decays, where D + s → φ(→ K + K − )π + and D 0 → K + π − . We apply the above 3D fit to B + → φφK + and B 0 → φφK 0 candidate events to determine the signal yield (and A CP in the first case). Figures 2 and 3 show M bc , ∆E, and O ′ NN projections of the fits. The fit results are listed in Table II. We find signal yields of 85.0 + 10.2 − 9.5 for B + → φφK + and 26.5 + 5.8 − 5.1 for B 0 → φφK 0 , and an A CP value of −0.02 ± 0.11 for the first case. We also apply the 3D fit to B + → φφK + candidate events with m φφ within the η c region to calculate the signal yield and A CP value. The corresponding M bc and ∆E projections are shown in Fig. 4, with the fit results listed in Table II. We obtain a signal yield of 73.2 + 9.0 − 8.3 and an A CP value of +0.12 ± 0.12 in the η c region. The signal significance is calculated as −2 log(L 0 /L max ), where L 0 and L max are the likelihood values with the signal yield fixed to zero and for the nominal fit, respectively. We include systematic uncertainties that impact only the signal yield into the likelihood curve via a Gaussian convolution before calculating the final significance.
where N 0 , N 1 , and N 2 are the yields obtained in SR, SB1, and SB2, respectively; n s , n a , and n b are the B → φφK yield in SR, B → φKKK yield in SB1, and B → KKKKK yield in SB2, respectively. Lastly, r s1 and r s2 are the ratios of B → φφK yields in SB1 and SB2 to that in SR; r a0 and r a2 are the ratios of B → φKKK yields in SR and SB2 to that in SB1; and r b0 and r b1 are the ratios of B → KKKKK yields in SR and SB1 to that in SB2. All these ratios are obtained from an MC study. We obtain the resonant B → φφK yield in SR (n s ) as 81.8 + 10.1 − 9.4 and 23.7 + 5.7 − 5.0 for the charged and neutral mode, respectively. These n s values are used in the branching fraction calculation of Eq. (5).
The background-subtracted distributions [18] of m φφ and m φK obtained for B ± → φφK ± below the η c threshold are shown in Fig. 6. These are broadly compatible with the predictions of a three-body phase space MC Events/(10 MeV) sample. In particular, we do not find any enhancement in the m φφ spectrum, including the 2.3 GeV/c 2 region [3] where a glueball and X(2350) candidates are predicted.
Systematic uncertainties in the branching fraction are listed in Table III. The uncertainties due to PDF shapes are estimated by varying all the fixed shape parameters by their errors. In particular, for fixed signal shape parameters, we vary the data-MC corrections by their uncertainties as determined using the control sample of TABLE II. Number of candidate events (n cand ), detection efficiency (ε), total and resonant signal yield (nsig), significance, branching fraction (B) and CP asymmetry (ACP ) obtained from a fit to data for B → φφK decays below and within the ηc region. Quoted uncertainties are statistical only, and significances defined in the text are given in terms of standard deviations.
Potential fit bias is checked by performing an ensemble test comprising 1000 pseudoexperiments, where signal is taken from the corresponding MC sample, and the PDF shapes are used to generate background events. We obtain a Gaussian normalized residual distribution of unit width, and add its mean and uncertainty in width in quadrature to calculate the systematic error. Uncertainty due to continuum suppression is obtained with the B + → D + s D 0 Events/(10 MeV)  control sample by comparing, between data and simulation, fit results obtained with and without the O NN requirement. A D * + → D 0 (K − π + )π + control sample is used to determine the systematic uncertainty due to the R K/π requirement. We use partially reconstructed D * + → D 0 (K 0 S π + π − )π + decays to assign the systematic uncertainty due to charged-track reconstruction (0.35% per track). The uncertainty due to K 0 S reconstruction is estimated from D 0 → K 0 S K 0 S decays [19]. We estimate the uncertainty due to efficiency variation across the Dalitz plot by weighting phase-space-generated signal MC events according to the measured distribution in data (Fig. 6) and taking the difference between the weighted and nominal efficiency. The total systematic uncertainty is obtained by adding all the above contributions in quadrature. We consider two possible sources of systematic uncertainties contributing to A CP , as listed in Table IV. The first is due to the intrinsic detector bias on charged kaon detection and is estimated using D + s → φπ + and D 0 → K − π + decays [20]. The second arises due to the potential variation of the PDF shapes. We calculate its contribution by following a procedure similar to that used in estimating the PDF shape uncertainties in the branching fractions.
In summary, we have measured the branching fractions and CP -violation asymmetries in B → φφK decays based on the full Υ (4S) data sample of 772 × 10 6 BB events collected by the Belle detector at the KEKB asymmetricenergy e + e − collider. We obtain the branching fraction and CP asymmetry for B ± → φφK ± below the η c threshold (m φφ < 2.85 GeV/c 2 ) as and −0.02 ± 0.11 ± 0.01, respectively. We also report the CP -violation asymmetry for B ± → φφK ± in the η c region (m φφ ∈ [2.94, 3.02] GeV/c 2 ) to be +0.12 ± 0.12 ± 0.01, consistent with no CP violation. The obtained value of the branching fraction of B ± → φφK ± decay is consistent with and supersedes our previous result [5]. The measured branching fraction for B 0 → φφK 0 below the η c threshold is