---
abstract: 'We present branching fraction and $CP$ asymmetry measurements as well as angular studies of $B \to \phi \phi K$ decays using $464\times 10^6$  events collected by the  experiment. The branching fractions are measured in the $\phi \phi$ invariant mass range below the $\eta_c$ resonance ($m_{\phi\phi}<2.85$ GeV). We find ${\cal B}(B^+ \to \phi \phi K^+) = (5.6 \pm 0.5 \pm 0.3)\times 10^{-6}$ and ${\cal B}(B^0 \to \phi \phi K^0) = (4.5 \pm 0.8 \pm 0.3)\times 10^{-6}$, where the first uncertaintiy is statistical and the second systematic. The measured direct $CP$ asymmetries for the $B^\pm$ decays are $A_{CP} = -0.10 \pm 0.08 \pm 0.02$ below the $\eta_c$ threshold ($m_{\phi\phi}<2.85$ GeV) and $A_{CP} = 0.09 \pm 0.10 \pm 0.02$ in the $\eta_c$ resonance region ($m_{\phi\phi}$ in \[2.94,3.02\] GeV). Angular distributions are consistent with $J^P = 0^-$ in the $\eta_c$ resonance region and favor $J^P = 0^+$ below the $\eta_c$ resonance.'
title: ' [ **Measurements of branching fractions and $CP$ asymmetries and studies of angular distributions for $B \to \phi \phi K$ decays**]{} '
---

BABAR-PUB-11/003; SLAC-PUB-14473\

The violation of $CP$ symmetry is a well-known requirement for the matter-antimatter imbalance of the universe [@sakarov]. The  [@babarnim] and Belle [@bellenim] experiments at the high-luminosity $B$ factories, PEP-II [@pepii] and KEKB [@kekb], have made numerous $CP$ asymmetry measurements using datasets two orders of magnitude larger than their predecessors. All of these measurements are consistent with a single source of $CP$ violation – the complex phase within the CKM quark mixing matrix of the Standard Model [@km]. However, with the small amount of $CP$ violation from the CKM matrix, it is difficult to explain the matter-antimatter asymmetry of the universe [@matter-antimatter-asym]. This motivates searches for new sources of $CP$ violation.

A method to search for new sources of $CP$-violating phases is to measure $CP$ asymmetries in hadron decays that are forbidden at the tree level [@grossman]. Since the leading decay amplitude is a one-loop process, contributions within the loop from virtual non-Standard-Model particles cannot be excluded. The quark interactions with the non-Standard-Model particles can introduce new $CP$ violating phases in the decay amplitude, which can lead to observable non-zero $CP$ asymmetries. Decays of $B$ mesons with a $b\to s \bar s s$ transition have been extensively studied for this reason.

The three body decay is a one-loop “penguin” $b\to s \bar s s$ transition. This final state can also occur through the tree-level decay , followed by , where the $B$ decay is a $b\to c\bar c s$ transition. If the $\phi\phi$ invariant mass in the three-body decay is close to the $\eta_c$ resonance, the tree and penguin amplitudes may interfere. Within the Standard Model, the relative weak phase between these amplitudes is $\arg(V_{tb}V_{ts}^*/V_{cb}V_{cs}^*) \approx 0$, so no $CP$ violation is expected from the interference. However, new physics contributions to the penguin loop in the decay could introduce a non-zero relative $CP$ violating phase, which may then produce a significant direct $CP$ asymmetry [@hazumi]. Measurement of a significant, non-zero direct $CP$ asymmetry would be an unambiguous sign of new physics. A previous measurement of the direct $CP$ asymmetry [@belle-ppk] was consistent with zero, but was also limited by a large statistical uncertainty. The $B^+$ and $B^0$ branching fractions have been previously measured [@belle-ppk][@old-babar-ppk] to be a few times $10^{-6}$. Theoretical predictions of the branching fractions are of the same order [@fajfer][@chen].

Dataset and Detector description
================================

We present measurements of the and branching fractions [@chargeconj] and direct $CP$ asymmetry $A_{CP} \equiv [N(B^-)-N(B^+)]/[N(B^-)+N(B^+)]$ as well as studies of angular distributions performed using $464 \times 10^6$ pairs collected by the  experiment at the SLAC National Accelerator Laboratory. The direct $CP$ asymmetry is measured both below and within the $\eta_c$ resonance region of the $\phi\phi$ invariant mass with these regions defined as $m_{\phi\phi}<2.85$ GeV and $m_{\phi\phi}$ within \[2.94, 3.02\] GeV, respectively [@c=1]. The branching fractions are measured in the region below the $\eta_c$ resonance ($m_{\phi\phi}<2.85$ GeV).

The  detector is described in detail elsewhere [@babarnim]. What follows is a brief overview of the main features of the detector. The detector has a roughly cylindrical geometry, with the axis along the beam direction. The trajectories, momenta, and production vertices of charged particles are reconstructed from position measurements made by a silicon vertex tracker (SVT) and a 40-layer drift chamber (DCH). The SVT consists of 5 layers of double-sided silicon strip detectors which provide precision position measurements close to the beam interaction region. Both the SVT and DCH measure the specific energy loss ($dE/dx$) along the charged particle trajectory, which is used to infer the particle mass from the velocity dependence of the energy loss and the momentum measurement. The tracking system is inside a uniform 1.5 T magnetic field provided by a superconducting solenoid. Outside the tracking system, an array of quartz bars coupled with an array of phototubes (DIRC) detects the Cherenkov light produced when a charged particle travels through the quartz bars. The measured Cherenkov angle is used to infer the particle mass from the velocity dependence of the Cherenkov angle and the measured momentum. The energies of photons and electrons are determined from the measured light produced in electromagnetic showers inside a CsI crystal calorimeter (EMC). Gaps in the iron of the magnet flux return are instrumented with resistive plate chambers and limited streamer tubes, which are used to identify muons and neutral hadrons (IFR).

We use Monte Carlo (MC) samples to determine the signal selection efficiency. The MC events are generated with EvtGen [@evtgen] and simulated using <span style="font-variant:small-caps;">Geant4</span> [@geant4].

Event selection
===============

We select events containing multiple hadrons by requiring at least three charged tracks in the event and the ratio of the second to zeroth Fox-Wolfram[@fox-wolfram] moments $R_2$ to be less than 0.98.

Charged kaon candidates are required to pass a selection based on a likelihood ratio which uses the SVT and DCH $dE/dx$ and the DIRC Cherenkov angle measurements as inputs to the likelihood. The ratio is defined as $R_{hh'} \equiv {\cal L}_h / ( {\cal L}_h + {\cal L}_{h'})$, where $h$ and $h'$ are $K$, $\pi$, or $p$. The minimum kaon selection criterion is $R_{K\pi}>0.2$ or $R_{p\pi}>0.2$. This selection has an efficiency greater than 98% for kaons and a pion efficiency of less than 15% below a lab momentum of 2.5 GeV. Candidate $\phi \to K^+K^-$ decays are constructed from oppositely-charged kaon candidates with an invariant mass in the range of 0.987 to 1.2 GeV. At least one of the kaons in each $\phi \to K^+K^-$ candidate must also satisfy the more stringent criteria of $R_{K\pi}>0.8176$ and $R_{Kp}>0.018$, which has an efficiency greater than 90% for kaons and a pion efficiency of less than 3% below a lab momentum of 2.5 GeV.

Candidate $K^0_S \to \pi^+\pi^-$ decays are constructed from oppositely-charged pion candidates with an invariant mass in the range of 0.486 to 0.510 GeV. The pion tracks are fit to a common vertex. The $\chi^2$ probability of the vertex fit must be greater than 0.001. The typical experimental resolution on the measured $K^0_S$ flight length in the plane transverse to the beam is around 0.2 mm or less. We require the transverse flight length to be at least 2 mm.

Candidate  decays are constructed from pairs of $\phi$ candidates that do not share any daughters and either a $K^\pm$ or a $K^0_S$ candidate. The $\phi$ and $K^\pm$ candidates are constrained to a common vertex. We reject combinatoric background by requiring the $B$ candidate to have kinematics consistent with $\Upsilon(4S) \to \BB$ using two standard variables: $m_{\rm ES}$ and $\Delta E$. The energy-substituted $B$ mass is defined as $m_{\rm ES} \equiv \sqrt{ E_{\rm beam}^{*2} - p^{*2} }$, where $E^*_{\rm beam}$ and $p^{*}$ are the beam energy ($\sqrt{s}/2$) and the reconstructed $B$ momentum, both in the center-of-mass (CM) reference frame. The energy difference is defined as $\Delta E \equiv E^* - E_{\rm beam}^{*}$, where $E^*$ is the reconstructed $B$ energy in the CM frame. We require $m_{\rm ES}$ and $\Delta E$ to be within \[5.20, 5.29\] GeV and \[$-0.1, 0.1$\] GeV, respectively. The experimental resolution is about 2.7 MeV for $m_{\rm ES}$ and 15 MeV for $\Delta E$.

The $m_{\rm ES}$ interval includes a large “sideband” region below the area where the signal is concentrated near the $B$ mass. The $\Delta E$ interval also is wide enough to include sideband regions where the signal probability is very low. Including events in the sideband regions enables us to determine the probability density functions (PDFs) of the combinatoric background directly in the maximum likelihood (ML) fits of the data.

About 7% of events in signal Monte Carlo samples have more than one  candidate. If there are multiple $B^+\to \phi\phi K^+$ candidates in a single event, we select the $B^+$ candidate with the smallest mass $\chi_m^2$ defined as $\sum_i \frac{(m_i - m_{0})^2}{\sigma^2}$, where the sum $i$ is over the two $\phi$ candidates, $m_i$ ($m_{0}$) is the reconstructed (nominal) $\phi$ mass, and $\sigma$ is the RMS of the reconstructed $m_{\phi}$ distribution for properly reconstructed $\phi$ candidates. If there are more than one $B^+$ candidates that use the same two $\phi$ candidates, we choose the $B^+$ candidate with the highest quality $K^+$ identification for the $K^+$ from the $B^+$ decay. If the quality level of the $K^+$ identification is the same for these $B^+$ candidates, we choose the $B^+$ candidate with the highest vertex $\chi^2$ probability. For events with multiple $B^0\to \phi\phi K^0_S$, the sum for $\chi_m^2$ includes the $K^0_S \to \pi^+\pi^-$, and we choose the $B^0$ with the smallest $\chi_m^2$. For both the $B^+$ and $B^0$ decay modes, the probability that the the algorithms described above choose the correct candidate is about 87%.

The reconstruction and selection efficiencies for events with $m_{\phi\phi} < 2.85$ GeV are determined from Monte Carlo samples to be 28.0% and 22.5% for the $B^+$ and $B^0$ modes, respectively.

We use control samples of $B \to D_s D$ decays where $D_s \to \phi \pi$, $D\to K\pi$, and $\phi \to K^+K^-$ to determine corrections to the $B\to \phi\phi K$ signal probability density function parameters determined from Monte Carlo samples in the maximum likelihood fits described below.

 Continuum Background 
----------------------

The events that pass the selection above with at least one $B$ candidate are primarily background events from the continuum ($e^+e^- \to q\bar q$ with $q = u,d,s,c$). We reduce this background by using a Fisher discriminant (${\cal F}$), which is the linear combination of seven variables and is optimized for maximum separation power of signal and the continuum background. The seven variables are listed below. These variables are commonly used by the  experiment in analyses of charmless $B$ decays, where the primary background is from continuum events. They take advantage of aspects of the production distributions and event topologies of versus continuum $q\bar q$ production events.\

-   [$|\Delta t/\sigma_{\Delta t}|$:]{} the absolute value of the reconstructed proper time difference between the two $B$ decays divided by its uncertainty [@dtandft].

-   [$|FT|$:]{} the absolute value of the standard  flavor tagging neural network output [@dtandft].

-   [$|\cos\theta^*_{\rm th}|$:]{} the absolute value of the cosine of the angle between the $B$ candidate thrust axis and the thrust axis of the rest of the event computed in the CM frame. The thrust axis is the direction that maximizes the scalar sum of the projection of the track momenta on that direction.

-   [$|\cos\theta^*_{B{\rm thr}}|$:]{} the absolute value of the cosine between the thrust axis of the $B$ candidate and the beam axis in the CM frame. Signal events have a uniform distribution in this variable, while continuum background follows a $|1+\cos^2\theta|$ distribution, where $\theta$ is the angle between the thrust direction and the beam axis.

-   [$|\cos\theta^*_B|$:]{} the absolute value of the cosine of the angle between the $B$ direction and the beam axis in the CM frame. The angular distribution of the signal follows a $\sin^2 \theta^*_B$ distribution, while the continuum background is uniformly distributed.

-   [$L_0$ and $L_2$]{}: The zeroth and second angular moments of the momentum flow of the rest of the event about the $B$ thrust axis, defined as $L_j \equiv \sum_i p_i |\cos \theta_i|^j$, where the angle $\theta_i$ is the angle between track $i$ and the $B$ thrust axis and the sum excludes the daughters of the $B$ candidate. The calculations are done in the CM frame.

Distributions of ${\cal F}$ for signal and continuum MC samples are shown in Fig. \[fig:fisher\]. The Fisher discriminant ${\cal F}$ is used as one of several variables in the maximum likelihood fits described below.

![ Distributions of the Fisher discriminant ${\cal F}$ for $B\to \phi \phi K$ signal (solid black) and $e^+e^- \to q\bar q$ continuum (dashed red) Monte Carlo simulation. []{data-label="fig:fisher"}](fisher1.eps){width="\linewidth"}

Peaking Backgrounds
-------------------

The ultimate detected state of our $B\to \phi\phi K$ signal decay is five kaons. In addition to the $\phi$ resonance, there may be contributions to each $K^+K^-$ pair either from other intermediate $K^+K^-$ resonances, such as the $f_0(980)$, or from non-resonant $K^+K^-$ contributions. We use the $K^+K^-$ mass sidebands for each $\phi$ candidate to determine the amount of $B$ mesons that decay to the detected five-kaon state (which we denote $B\to 5K$) that are not coming from $B\to \phi\phi K$. The specific $B$ decays that we consider as backgrounds are $B\to \phi K^+K^-K$, $B\to K^+K^-K^+K^-K$, $B\to f_0 \phi K$, and $B\to f_0 K^+K^-K$. The branching fractions for these decays are currently unknown. We call these $B$ decays “peaking backgrounds” because properly reconstructed $B$ candidates are indistinguishable from our $B\to \phi\phi K$ signal in the $m_{\rm ES}$, $\Delta E$, and ${\cal F}$ variables.

We perform unbinned extended maximum likelihood fits to determine the signal and combinatoric background yields and, in some cases, the charge asymmetry. All of the fits use the product of one-dimensional PDFs of $m_{\rm ES}$, $\Delta E$, and ${\cal F}$ in the likelihood. For the $B\to \phi\phi K$ branching fraction measurements, we also include PDFs for the invariant mass of each $\phi \to K^+K^-$ candidate ($m_{\phi 1}$ and $m_{\phi 2}$).

As a first step, we divide the $m_{\phi 1}$ vs. $m_{\phi 2}$ plane [@phi-assignment] in the range of 0.987 to 1.200 GeV into five mutually exclusive zones. We fit for the $B \to 5K$ yield in each zone using only $m_{\rm ES}$, $\Delta E$, and ${\cal F}$ in the likelihood. The zones are based on various combinations of the $\phi$ signal and sideband regions, which are defined as: [Low-SB]{} \[0.987,1.000\] GeV, [phi-signal]{} \[1.00,1.04\] GeV, and [High-SB]{} \[1.04,1.20\] GeV. Each of the five zones is chosen so that either the $B\to \phi\phi K$ signal or one of the four peaking $B$ backgrounds is concentrated in the region. We compute the number of peaking background events within the $m_{\phi}$ range used for the branching fraction fit by using the results of the five zone fits as described below.

Figure \[fig:mphi1vsmphi2\] shows the distribution of events in the $m_{\phi 2}$ vs. $m_{\phi 1}$ plane for the selected $B^+\to 5K$ candidates in the data. To enhance the $B^+\to 5K$ signal for the figure, we have required $m_{\rm ES}>5.27$ GeV, $|\Delta E|<0.040$ GeV, and ${\cal F}<0.0$. The inset of the figure shows the definition of the five zones. A concentration of events in the [phi-signal]{} region for both $\phi$ candidates (zone 1) is clearly evident. The region defined as [phi-signal]{} combined with [High-SB]{} for either $\phi$ candidate (zone 2) contains the largest fraction of the $B\to \phi K^+K^-K$ mode, although the $B\to \phi\phi K$ signal also populates this region due to cases where one $\phi$ is mis-reconstructed. The zone where the invariant mass of both $\phi$ candidates is in the [High-SB]{} region (zone 3) contains the largest concentration of the non-resonant $B\to K^+K^-K^+K^-K$ mode. Zones 4 and 5 contain a large fraction of the $B\to f_0 \phi K$ and $B\to f_0 K^+K^- K$ modes, respectively, and very small fractions of the other three modes.

![ Data distribution of $B^+ \to 5K$ events in the $m_{\phi 2}$ vs. $m_{\phi 1}$ plane for events with $m_{\phi\phi} < 2.85$ GeV. To enhance the $B\to 5K$ signal for the figure, we have required $m_{\rm ES}>5.27$ GeV, $|\Delta E|<0.040$ GeV, and ${\cal F}<0.0$. The efficiency of these additional requirements, relative to the nominal selection, is about 70% for the signal. The inset shows the definition of the five zones. []{data-label="fig:mphi1vsmphi2"}](mphi1vsmphi2sp.eps){width="\linewidth"}

Monte Carlo samples for the five $B$ decay modes (signal plus four peaking background modes) are used to determine the fraction of events in each zone ($i$) for each decay mode ($j$), which we denote with the matrix $f_{ij}$. The total $B\to 5K$ yield ($n_i$) is determined for each zone $i$ using five separate maximum likelihood fits of the data. The yield for each $B$ decay mode ($N_j$) and the amount of each mode $j$ in zone $i$ ($n_{ij}$) can be determined from $$N_j = \sum_i f^{-1}_{ij} n_i \ \ \ \ {\rm and} \ \ \ n_{ij} = f_{ij} N_j .$$ Zone 1 corresponds to the $m_{\phi}$ range used in the branching fraction maximum likelihood fit.

 Maximum Likelihood Fits 
-------------------------

The extended maximum likelihood fits in the five zones determine the $B\to 5K$ signal and combinatoric background yields in each zone. The $B\to 5K$ signal is split into properly reconstructed and misreconstructed (“self-crossfeed”) components, with the self-crossfeed fraction fixed. The self-crossfeed component is defined as events where a true $B\to 5K$ decay is present in the event, but one or more tracks used in the reconstructed $B$ are either from the other $B$ in the event or not real. In zone 1, the self-crossfeed fraction for $B\to \phi\phi K$ decays is around 7%. The properly reconstructed $B\to 5K$ signal component is described by the following PDFs: a Crystal Ball function [@CB-function] for $m_{\rm ES}$, the sum of three Gaussians for $\Delta E$, and the sum of a bifurcated Gaussian and a Gaussian for ${\cal F}$. The Crystal Ball function is a Gaussian modified to have an extended power-law tail on the low side. The $B\to 5K$ signal PDF parameters are determined from MC samples with corrections to the $m_{\rm ES}$ and $\Delta E$ core mean and width parameters from the $B\to D_s D$ control samples. The mean corrections are $0.04 \pm 0.11$ MeV and $-3.5 \pm 0.8$ MeV for $m_{\rm ES}$ and $\Delta E$, respectively. The width scale factors are $1.10 \pm 0.04$ and $1.04 \pm 0.05$ for $m_{\rm ES}$ and $\Delta E$, respectively. The combinatoric background is described by the following PDFs: an empirical threshold function [@Argus-function] for $m_{\rm ES}$, a first-order polynomial for $\Delta E$, and the sum of two Gaussians for ${\cal F}$. Most of the combinatoric background PDF shape parameters are determined in the fits.

The results of the five zone fits for the $B^+$ and $B^0$ modes are given in Tables \[tab:b+zone-fits\] and \[tab:b0zone-fits\], respectively, in the appendix. The $B\to \phi\phi K$ signal is observed in both the $B^+$ and $B^0$ samples. The $B^0 \to \phi\phi K^0$ decay has not been observed previously. The $B$ yield in zone 2 for the $B^+$ mode is significant, but about half of this is due to misreconstructed $B\to \phi\phi K$ signal. The computed $B^+\to \phi K^+K^-K^+$ and $B^+\to K^+K^-K^+K^-K^+$ yields are positive, but the significance is less than two standard deviations. There is no evidence of either $B\to f_0\phi K$ or $B \to f_0 K^+K^-K$. The branching fraction maximum likelihood fits use the $m_{\phi}$ range that corresponds to zone 1. We fix the yield of each of the four peaking background modes to the zone 1 value in Table \[tab:b+zone-fits\] or \[tab:b0zone-fits\] for the branching fraction fit described below.

Branching Fraction Analysis 
============================

The maximum likelihood fit used to measure the $B\to \phi \phi K$ yield below the $\eta_c$ resonance for the branching fraction measurement restricts the event selection with $m_{\phi\phi} < 2.85$ GeV and $m_{\phi}$ within \[1.00,1.04\] GeV, which corresponds to zone 1 in the peaking background discussion above. The fit components are $B\to \phi \phi K$ signal, combinatoric background, and the four peaking backgrounds.

In addition to $m_{\rm ES}$, $\Delta E$, and ${\cal F}$, PDFs for $m_{\phi 1}$ and $m_{\phi 2}$ are included in the likelihood function. For each fit component, each $\phi$ candidate has a PDF that is the sum of a properly reconstructed $\phi \to K^+K^-$ decay, given by a relativistic Breit-Wigner function, and a misreconstructed $\phi$, described by a first-order polynomial. The $m_{\phi 1}$ and $m_{\phi 2}$ PDFs are combined in a way that is symmetric under $1 \leftrightarrow 2$ exchange and takes into account the fractions of events where both $\phi$ candidates are properly reconstructed, one $\phi$ is misreconstructed, and both $\phi$ candidates are misreconstructed.

In addition to the signal and combinatoric background yields, the charge asymmetry for the signal and combinatoric background components and most of the combinatoric background PDF parameters are determined in the fit.

The results of the $B^+$ and $B^0$ fits are shown in Figs. \[fig:bplus-fit\] and \[fig:b0-fit\], respectively. To reduce the combinatoric background in each distribution shown in the figures, a requirement is made on a likelihood ratio, which is based on all the fit variables except the one plotted.

![ Results of fitting the $B^+ \to \phi\phi K^+$ sample for $m_{\phi\phi}<2.85$ GeV. The dashed red curve is the sum of the combinatoric and peaking background components. The solid blue curve is for all components. A requirement on a likelihood ratio based on all fit variables except the one plotted is made to reject most of the background. The likelihood ratio requirements are about 84% efficient for the signal. []{data-label="fig:bplus-fit"}](bplus-fit.eps){width="\linewidth"}

![ Results of fitting the $B^0 \to \phi\phi K^0_S$ sample for $m_{\phi\phi}<2.85$ GeV. The dashed red curve is the sum of the combinatoric and peaking background components. The solid blue curve is for all components. A requirement on a likelihood ratio based on all fit variables except the one plotted is made to reject most of the background. The likelihood ratio requirements are about 84% efficient for the signal. []{data-label="fig:b0-fit"}](b0-fit.eps){width="\linewidth"}

The fitted charge asymmetry ${\cal A}_{CP}$ for the background component is $0.02 \pm 0.03$. The charge asymmetry for the signal component is $-0.10 \pm 0.08$. The fitted yields of $B^+$ and $B^0$ signal candidates with $m_{\phi\phi}<2.85$ GeV are $178 \pm 15$ events and $40 \pm 7$ events, respectively, where the uncertainties are statistical only.

Systematic Uncertainties
------------------------

Table \[tab:bf-fit-syst\] summarizes the systematic uncertainties on the $B\to \phi\phi K$ branching fractions in the $m_{\phi\phi}<2.85$ GeV region. The systematics are divided into additive uncertainties that affect the $B$ yield measurement and multiplicative uncertainties in the branching fraction calculation.

The uncertainties from the corrections applied to the PDF parameters such as the $m_{\rm ES}$ and $\Delta E$ core mean and width for the signal component, which are derived from data control samples, are listed under “ML Fit Yield”. The signal Fisher and $m_\phi$ core Gaussian mean and width parameters are not corrected in the fit, because data control sample measurements are consistent with the Monte Carlo. However, we did vary the signal Fisher and $m_\phi$ core Gaussian mean and width parameters by the statistical uncertainty of the data control sample measurements. These variations are also included under “ML Fit Yield”. The fit bias systematic is taken to be half of the bias correction added in quadrature with the statistical uncertainty on the bias. We vary the fixed peaking background yields by their statistical uncertainties (see Tables \[tab:b+zone-fits\] and \[tab:b0zone-fits\]) and by varying the fractions $f_{ij}$. The fixed self-crossfeed fraction for the signal component was varied by $\pm 2$%. Adding the individual uncertainties in quadrature, the total additive systematic uncertainties on the $B^+$ and $B^0$ signal yields are 6.2 and 1.8 events, respectively.

The uncertainty on the track reconstruction efficiency is $\pm 0.23$% per track, which is taken to be fully correlated for the charged kaons. The $K^0_S$ reconstruction efficiency has an uncertainty of 1.5%. The $\phi\to K^+K^-$ and $K^0_S\to \pi^+\pi^-$ branching fractions are taken from the PDG [@pdg] and are varied by their one standard deviation uncertainties. The systematic uncertainty on the $K^\pm$ identification criteria was estimated by comparing the ratio of the $B$ yield with the nominal selection to the $B$ yield requiring all $K^\pm$ to pass the tighter selection in the data and the MC samples. This gives an uncertainty of 3% for the $B^+$ mode and 2% for the $B^0$ mode. Adding the individual uncertainties in quadrature, the overall multiplicative systematic uncertainties are 3.6% for the $B^+$ mode and 3.2% for the $B^0$ mode.

[ccc]{} Quantity & $B^+ \rightarrow \phi\phi K^+$ & $B^0 \rightarrow \phi\phi K^{0}$\
Fit Stat. Uncertainty (events) & 15.1 & 7.0\
\
ML Fit Yield & 3.3 & 1.0\
ML Fit Bias & 1.6 & 0.2\
Peaking BG, region fits & 3.5 & 1.2\
Peaking BG, $f_{ij}$ values & 3.0 & 0.8\
Self-Crossfeed Fraction & 1.8 & 0.4\
Total Additive Syst. (events)& 6.2 & 1.8\
\
Tracking Efficiency & 1.2 & 1.0\
$K_s^0$ Reconstruction Efficiency & - & 1.5\
Number & 1.1 & 1.1\
$\cal B$ ($\phi \rightarrow K^+ K^-$) & 1.2 & 1.2\
$\cal B$ ($K_s^0 \rightarrow \pi^+ \pi^-$) & - & 0.1\
MC Statistics & 0.1 & 0.1\
$m_{\phi\phi}$ Cut Efficiency & 0.2 & 0.3\
$K^\pm$ Identification & 3.0 & 2.0\
Total Multiplicative Syst. (%) & 3.6 & 3.2\
Total Systematic \[$\cal B$\] ($\times 10^{-6}$) & 0.3 & 0.3\
Statistical \[$\cal B$\] ($\times 10^{-6}$) & 0.5 & 0.8\

The signal charge asymmetry has been corrected for a bias due to differences in the $K^+$ and $K^-$ efficiencies by adding $+0.010\pm0.005$ to the asymmetry. The overall 2% systematic uncertainty takes into account uncertainties on the charge dependence of the tracking efficiency, material interaction cross section for kaons, and particle identification.

Branching Fraction Results
--------------------------

Table \[tab:bf-fit-results\] summarizes the $B \to \phi\phi K$ branching fraction results for $m_{\phi\phi} < 2.85$ GeV. We find $$\begin{aligned}
     {\cal B}(B^+ \to \phi\phi K^+) & = & (5.6 \pm 0.5 \pm 0.3) \times 10^{-6} \\
     {\cal B}(B^0 \to \phi\phi K^0) & = & (4.5 \pm 0.8 \pm 0.3) \times 10^{-6},
   \end{aligned}$$ where the first uncertainty is statistical and the second systematic. These results are consistent with and supersede the previous measurements [@old-babar-ppk] by the  Collaboration. The Belle collaboration measurements [@belle-ppk] are lower, though they are statistically compatible. Our branching fraction measurements are higher than the theoretical predictions of [@fajfer] and [@chen].

                             $B^+ \rightarrow \phi\phi K^+$   $B^0 \rightarrow \phi\phi K^{0}$
  ------------------------- -------------------------------- ----------------------------------
        Events to fit                     1535                              293
      Fit signal yield                178 $\pm$ 15                       40 $\pm$ 7
    ML-fit bias (events)             3.0 $\pm$ 0.5                     0.0 $\pm$ 0.2
      MC efficiency (%)                   28.0                              22.5
     $\Pi{\cal B}_i$(%)                   24.2                              8.4
     Stat. significance                    24                                11
                                                             
     $\cal B$($10^{-6}$)        5.6 $\pm$ 0.5 $\pm$ 0.3           4.5 $\pm$ 0.8 $\pm$ 0.3
                                                             
     Signal [$A_{CP}$]{}       $-0.10 \pm 0.08 \pm 0.02$                     -
   Comb. Bkg. [$A_{CP}$]{}          0.02 $\pm$ 0.03                          -

  :  Branching fraction and charge asymmetry results for $B\to \phi\phi K$ in the region $m_{\phi\phi} < 2.85$ GeV. The statistical significance is given by $\sqrt{2 \ln ({\cal L}_{\rm max}/{\cal L}_0)}$, where ${\cal L}_{\rm max}$ is the maximum likelihood and ${\cal L}_0$ is the likelihood for the hypothesis of no $\phi\phi K$ signal. The significance does not include systematic uncertainties. []{data-label="tab:bf-fit-results"}

CP asymmetry in $\eta_c$ resonance region
=========================================

As was mentioned in the introduction, a significant non-zero direct $CP$ asymmetry in the $\eta_c$ resonance region of $m_{\phi\phi}$ would be a clear sign of physics beyond the Standard Model. For this measurement, we use the simpler likelihood, based on $m_{\rm ES}$, $\Delta E$, and ${\cal F}$. Figure \[fig:mphiphi\] shows the fitted $B^+\to \phi\phi K^+$ yield as a function of $m_{\phi\phi}$. The $\eta_c$ resonance is clearly visible. Narrow bins around the $\chi_{c0}$ and $\chi_{c2}$ resonances do not show a significant excess above the broad non-resonant component.

![ Fitted $B^+\to \phi\phi K^+$ yield as a function of $m_{\phi\phi}$. Each point shows the results of a maximum likelihood fit of the events in that bin. The inset is the same data with an expanded vertical range to show the shape of the non-resonant component more clearly. The yield has been divided by the bin width and scaled by 0.027 GeV, which is the bin width of the three bins in the $\eta_c$ resonance region (\[2.94,3.02\] GeV and dashed vertical lines in the inset). The two narrow bins above the $\eta_c$ are centered on the $\chi_{c0}$ (bin range \[3.400,3.430\] GeV) and the $\chi_{c2}$ (bin range \[3.552,3.560\] GeV). []{data-label="fig:mphiphi"}](prd-mphiphi.eps){width="\linewidth"}

The results of fitting the events in the $m_{\phi\phi}$ range of \[2.94,3.02\] GeV are given in Table \[tab:etac-region-fit-results\]. For the $CP$ asymmetry, we find $$A_{CP}(m_{\phi\phi}\ {\rm in} \ [2.94,3.02] \ {\rm GeV}) = -0.09 \pm 0.10 \pm 0.02 ,$$ where the first uncertainty is statistical and the second uncertainty is systematic. The value above includes the same 1% bias correction and has the same 2% overall systematic uncertainty as the signal charge asymmetry below the $\eta_c$ resonance as described above.

The fit yields $100 \pm 10$ signal candidates. Using ${\cal B}(B^+\to \eta_c K^+) = (9.1 \pm 1.3) \times 10^{-4}$ and ${\cal B}(\eta_c \to \phi\phi) = (2.7 \pm 0.9) \times 10^{-3}$ from the PDG [@pdg], a $B^+\to \phi\phi K^+$; $\phi\to K^+K^-$ reconstruction efficiency of 29% in the $\eta_c$ resonance region, and an efficiency of 78% for the $m_{\phi\phi}$ window of \[2.94,3.02\] GeV for the $\eta_c$ resonance, we would expect $62 \pm 22$ signal events, ignoring the non-resonant $B^+\to \phi\phi K^+$ contribution and any interference between the resonant $\eta_c$ and non-resonant amplitudes. We do not use our $B^+$ event yield to measure ${\cal B}(B^+\to \eta_c K^+) \times {\cal B}(\eta_c \to \phi\phi)$ due to the potentially large interference effects between the resonant and non-resonant $\phi\phi$ amplitudes which we can not easily quantify.

The $A_{CP}$ may integrate to zero, even if there is a contributing non-Standard-Model amplitude with a non-zero $CP$ violating phase. However, in this case the phase variation of the $\eta_c$ resonance amplitude could give non-zero $A_{CP}$ values with opposite signs above and below the peak of the resonance. We have performed the measurement in two ranges, splitting the $\eta_c$ region into two regions (above and below the peak of the resonance). The results are $$\begin{aligned}
      A_{CP}(m_{\phi\phi}\ {\rm in} \ [2.94,2.98] \ {\rm GeV}) & = & -0.10 \pm 0.15 \pm 0.02 \\
      A_{CP}(m_{\phi\phi}\ {\rm in} \ [2.98,3.02] \ {\rm GeV}) & = & -0.08 \pm 0.14 \pm 0.02,
   \end{aligned}$$ both of which are consistent with zero, as expected in the Standard Model.

   ML fit quantity/Analysis    $B^+ \rightarrow \phi\phi K^+$
  --------------------------- --------------------------------
         Events to fit                      181
       Fit signal yield                 100 $\pm$ 10
       MC efficiency (%)                    29.2
   Corr. Signal [$A_{CP}$]{}     $-0.09 \pm 0.10 \pm 0.02$
    Comb. Bkg. [$A_{CP}$]{}           $-0.06 \pm 0.11$

  : Fit results for $B^+\to \phi\phi K^+$ within $\eta_c$ resonance region ($m_{\phi\phi}$ within \[2.94,3.02\] GeV). The signal charge asymmetry [$A_{CP}$]{} has been corrected by adding $+0.010 \pm 0.005$ to the fitted asymmetry. []{data-label="tab:etac-region-fit-results"}

Angular Studies
===============

We use the angular variables that describe the $B^+\to \phi\phi K^+$ decay to investigate the spin components of the $\phi\phi$ system below and within the $\eta_c$ resonance. The angles are defined as follows.

-   [$\theta_i$]{}, $(i=1,2)$ : The $\theta_i$ angle is the angle between the momentum of the $K^+$ coming from the decay of $\phi_i$ in the $\phi_i$ rest frame with respect to the boost direction from the $\phi\phi$ rest frame to the $\phi_i$ rest frame.

-   [$\chi$]{}: The $\chi$ angle is the dihedral angle between the $\phi_1$ and $\phi_2$ decay planes in the $\phi\phi$ rest frame.

-   [$\theta_{\phi\phi}$]{}: The $\theta_{\phi\phi}$ angle is the angle between one of the $\phi$ mesons in $\phi\phi$ rest frame with respect to the boost direction from the $B^+$ rest frame to the $\phi\phi$ rest frame.

We project the $J^P = 0^-$ component by making a histogram of $m_{\phi\phi}$ weighting each event by $$\begin{split}
           P_2(\cos \theta_1) \ {\rm Re} \left[ Y^2_2( \theta_2, \chi) \right] \ = \\
                   \frac{25}{4} \left\{ 3 \cos^2 \theta_1  - 1 \right\}
                            \, \sin^2 \theta_2  \, \cos  2 \chi,
   \end{split}$$ where $P_2$ is a second-degree Legendre polynomial and $Y^2_2$ is a spherical harmonic with $\ell=2$ and $m=2$. In each bin, the $J^P = 0^-$ component yield is projected out, while the combinatoric background averages to zero. To do this, we select events in a signal region defined by: $m_{\rm ES} > 5.27$ GeV, $|\Delta E|< 40$ MeV, $m_\phi$ within \[1.01,1.03\] GeV, and ${\cal F}<0.5$. The efficiency of these requirements, relative to the selection used in the asymmetry measurement, is about 78% for signal events and 2.9% for combinatoric background. The combinatoric background that remains after this selection is shown using data events in the sideband region ($m_{\rm ES} <5.27$ GeV and $|\Delta E|<100$ MeV) scaled by 0.065, which is the signal-to-sideband ratio for the combinatoric background.

The results are shown in Fig. \[fig:angular-weight\]. The weighted yield in the $\eta_c$ region is consistent with all of the $B^+\to \phi\phi K^+$ events having $J^P = 0^-$. Just below the $\eta_c$ region, the weighted yield is consistent with zero. The excess in the bins near 2.2 GeV may be due to the $\eta(2225)$ seen in $J/\psi \to \gamma \phi\phi$ events at Mark III [@markIII] and BES [@BES].

Figure \[fig:angular-plots\] shows background-subtracted distributions of $\chi$, $\cos\theta_i$, and $|\cos\theta_{\phi\phi}|$ for the nominal event selection. The background subtraction is done with the technique described in reference [@splots]. Since there is no meaningful distinction between $\phi_1$ and $\phi_2$, we combine the $\cos\theta_1$ and $\cos\theta_2$ distributions into one plot of $\cos\theta$. The reconstruction and selection efficiency, determined from MC samples, is uniform in $\chi$ and $\cos\theta_1$, but not in $|\cos\theta_{\phi\phi}|$, so the $|\cos\theta_{\phi\phi}|$ distribution is efficiency corrected. For each distribution, we performed a simple least-$\chi^2$ fit to the distributions expected for both $J^P = 0^-$ and $J^P = 0^+$ for the $\phi\phi$ system.

![ Histograms (top) and weighted distributions (bottom) of $m_{\phi\phi}$ for the signal region (solid points) and data sideband selection (dashed with open diamonds) defined in the text. The sideband distributions have been normalized to the expected level of combinatoric background remaining after the signal region selection. Events in the bottom distribution were weighted by $P_2(\cos \theta_1) \ {\rm Re} \left[ Y^2_2( \theta_2, \chi) \right]$ which projects out the $J^P = 0^-$ component. The yield has been divided by the bin width and scaled by 0.211 GeV, which is the width of the bin covering the $\eta_c$ resonance (\[2.875,3.086\] GeV). []{data-label="fig:angular-weight"}](mphiphi-angularweight-a.eps "fig:"){width="0.99\linewidth"} ![ Histograms (top) and weighted distributions (bottom) of $m_{\phi\phi}$ for the signal region (solid points) and data sideband selection (dashed with open diamonds) defined in the text. The sideband distributions have been normalized to the expected level of combinatoric background remaining after the signal region selection. Events in the bottom distribution were weighted by $P_2(\cos \theta_1) \ {\rm Re} \left[ Y^2_2( \theta_2, \chi) \right]$ which projects out the $J^P = 0^-$ component. The yield has been divided by the bin width and scaled by 0.211 GeV, which is the width of the bin covering the $\eta_c$ resonance (\[2.875,3.086\] GeV). []{data-label="fig:angular-weight"}](mphiphi-angularweight-b.eps "fig:"){width="0.99\linewidth"}

![image](prd-chi-splot-etac-nostats.eps){width="0.32\linewidth"} ![image](prd-costheta-splot-etac-nostats.eps){width="0.32\linewidth"} ![image](prd-costhpp-etac-nostats.eps){width="0.32\linewidth"} ![image](prd-chi-splot-bec-nostats.eps){width="0.32\linewidth"} ![image](prd-costheta-splot-bec-nostats.eps){width="0.32\linewidth"} ![image](prd-costhpp-bec-nostats.eps){width="0.32\linewidth"}

   $m_{\phi\phi}$ (GeV)            Variable                 PDF        $\chi^2/N_{\rm dof}$     $\chi^2$ prob.
  ---------------------- ---------------------------- --------------- ---------------------- ---------------------
      $[2.94, 3.02]$                $\chi$             $\sin^2 \chi$        $9.51/9$                $0.39$
      $[2.94, 3.02]$                $\chi$                uniform            $60.3/9$         $1.2\times 10^{-9}$
         $<2.85$                    $\chi$             $\sin^2 \chi$         $41.6/9$         $3.9\times 10^{-6}$
         $<2.85$                    $\chi$                uniform            $18.5/9$               $0.030$
      $[2.94, 3.02]$            $\cos \theta$          $\sin^2 \chi$         $9.97/9$               $0.39$
      $[2.94, 3.02]$            $\cos \theta$             uniform            $60.5/9$         $1.1\times 10^{-9}$
         $<2.85$                $\cos \theta$          $\sin^2 \chi$         $32.9/9$         $1.7\times 10^{-4}$
         $<2.85$                $\cos \theta$             uniform            $25.8/9$         $2.2\times 10^{-3}$
      $[2.94, 3.02]$      $|\cos \theta_{\phi\phi}|$      uniform            $9.02/9$               $0.44$
         $<2.85$          $|\cos \theta_{\phi\phi}|$      uniform            $5.01/9$               $0.83$

  : Quality of the angular fits shown in Fig. \[fig:angular-plots\]. The first column is the $m_{\phi\phi}$ interval for the events in the fit. The last column is the $p$-value of the $\chi^2$ goodness-of-fit test for the hypothesis indicated in the third column.[]{data-label="tab:angular-fits"}

For a $J^P = 0^-$ state, we expect $\chi$ to have a $\sin^2 \chi$ distribution, while $\chi$ should be uniform for $J^P = 0^+$. The signal events in the $\eta_c$ resonance region are consistent with a $\sin^2 \chi$ distribution while the signal below the $\eta_c$ resonance is not. For a $J^P = 0^-$ state, the distributions of $\cos \theta_i$ are expected to have $\sin^2 \theta_i$ distributions, while a $J^P = 0^+$ state is expected to have uniform $\cos \theta_i$ distributions. The events in the $\eta_c$ resonance region are consistent with a $\sin^2 \theta_i$ distribution, while the events below the $\eta_c$ resonance show a deviation from a $\sin^2 \theta_i$ shape.

Finally, a spin-zero state should have a uniform $|\cos \theta_{\phi\phi} |$ distribution. The efficiency-corrected distributions shown in Fig. \[fig:angular-plots\], both within and below the $\eta_c$ resonance region, are consistent with a uniform $|\cos \theta_{\phi\phi} |$ distribution.

Summary and Conclusions
=======================

We have measured the branching fractions and charge asymmetries of $B\to \phi\phi K$ decays below the $\eta_c$ resonance in the $\phi\phi$ invariant mass ($m_{\phi\phi} < 2.85$ GeV). We observe both $B^+ \to \phi\phi K^+$ and $B^0 \to \phi\phi K^0_S$, each with a significance of greater than five standard deviations. The $B^0 \to \phi\phi K^0_S$ decay has not been observed previously. Our branching fraction measurements are higher than the theoretical predictions of [@fajfer] and [@chen].

We have measured the charge asymmetry for $B^+\to \phi\phi K^+$ in the $\eta_c$ resonance region, where a significant non-zero value would be an unambiguous indication of new physics. Our measurement is consistent with zero, which is the expectation of the Standard Model.

Finally, we have studied the angular distributions of $B^+ \to \phi\phi K^+$ decays below and within the $\eta_c$ resonance. We conclude from these studies that the non-resonant $B^+ \to \phi\phi K^+$ events below the $\eta_c$ resonance are, on average, more consistent with $J^P = 0^+$ than $J^P = 0^-$, while the distributions within the $\eta_c$ resonance region are all consistent with $J^P = 0^-$.

 Peaking background zone fit results
====================================

The results of the fits of the $B\to 5K$ yield in the five zones in the $m_{\phi 2}$ vs $m_{\phi 1}$ plane (shown in Fig. \[fig:mphi1vsmphi2\]) and the derived yields for each of the five $B$ decay modes are given in Tables \[tab:b+zone-fits\] and  \[tab:b0zone-fits\] below.

   Zone      $B^+$ yield                                                                         
  ------- ------------------ ------------------ --------------- --------------- ---------------- ---------------
     1     $188.4 \pm 16.0$   $183.2 \pm 17.0$   $6.3 \pm 3.9$   $0.8 \pm 0.8$    $0 \pm 3.0$     $0.3 \pm 1.5$
     2     $ 84.4 \pm 18.0$          36               39               9              -1.4              2
     3     $ 49.7 \pm 19.0$         3.8               18              26              -0.2             2.2
     4     $  1.0 \pm  2.0$         1.3               0.9             0.2             -1.7             0.3
     5     $  3.5 \pm  5.0$         0.2               0.8             1.3             -0.2             1.3
   $1-5$                        $225 \pm 21$      $65 \pm 40$     $38 \pm 38$    $-5.7 \pm 7.6$   $6.2 \pm 26$

   Zone     $B^0$ yield                                                                      
  ------- ---------------- ---------------- --------------- --------------- ---------------- ----------------
     1     $43.4 \pm 8.0$   $42.8 \pm 8.4$   $1.2 \pm 1.5$   $0.1 \pm 0.3$    $0 \pm 0.7$      $0 \pm 0.02$
     2     $15.1 \pm 7.0$         8                7              0.9             -0.4             -0.4
     3     $ 6.0 \pm 6.0$         1                3              2.5             -0.1             -0.4
     4      $ 0 \pm 0.5$         0.3              0.2              0              -0.5              0
     5      $ 0 \pm 0.8$          0               0.1             0.1              0               -0.2
   $1-5$                     $52 \pm 10$      $12 \pm 15$    $3.6 \pm 12$    $-1.6 \pm 1.8$   $-1.1 \pm 4.3$

[99]{} A.D. Sakharov, JETP Lett. [**5**]{}, 24 (1967). B. Aubert [*et al.*]{} \[The  Collaboration\], Nucl. Instrum. Meth. [**A479**]{}, 1 (2002). A. Abashian [*et al.*]{} \[The Belle Collaboration\], Nucl. Instrum. Meth. [**A479**]{}, 117 (2002). PEP-II Conceptual Design Report, SLAC-0418 (1993). S. Kurokawa and E. Kikutani, Nucl. Instrum. Meth. [**A499**]{}, 1 (2003). M. Kobayashi and T. Maskawa, Prog. Theor. Phys. [**49**]{}, 652 (1973). V. A. Rubakov and M. E. Shaposhnikov, Phys. Usp. [**39**]{}, 461 (1996). Y. Grossman and M. P. Worah, Phys. Lett. [**B 395**]{}, 241 (1997). M. Hazumi, Phys. Lett. [**B 583**]{}, 285 (2004). H. C. Huang [*et al.*]{} \[The Belle Collaboration\], Phys. Rev. Lett. [**91**]{}, 241802 (2003).; Y. T. Shen, K. F. Chen, P. Chang [*et al.*]{} \[The Belle Collaboration\], arXiv:0802.1547\[hep-ex\]. B. Aubert [*et al.*]{} \[The  Collaboration\], Phys. Rev. Lett. [**97**]{}, 261803 (2006). S. Fajfer, T. N. Pham, and A. Prapotnik, Phys. Rev. [**D69**]{}, 114020 (2004). C.-H. Chen and H.-n. Li, Phys. Rev. [**D70**]{}, 054006 (2004). Charge-conjugate states are implicitly included throughout the paper unless stated otherwise. Throughout the paper, necessary factors of $c$ are implied in the units for energy, mass, and momentum. D. Lange, Nucl. Instrum. Meth., A [**462**]{}, 152 (2001). S. Agostinelli [*et al.*]{} \[The <span style="font-variant:small-caps;">Geant4</span> Collaboration\], Nucl. Instrum. Meth., A [**506**]{}, 250 (2003). G. Fox and S. Wolfram, Phys. Rev. Lett., [**41**]{}, 1581 (1978). B. Aubert [*et al.*]{} \[The  Collaboration\], Phys. Rev. [**D79**]{}, 072009 (2009). In each event one of the $\phi$ candidates is randomly chosen to be $\phi_1$ with the other $\phi_2$. All aspects of the analysis, such as the PDFs used in the likelihood fits, are symmetric under the exchange of the $\phi_1$ and $\phi_2$ assignments (e.g. $m_{\phi 1} \leftrightarrow m_{\phi 2}$). J. E. Gaiser, Ph.D. thesis, Stanford University \[SLAC-R-255\] (1982). With $x\equiv m_{\rm ES}/E_{\rm beam}^{*}$ and $\xi$ a parameter to be fitted, $f(x) \propto x \sqrt{1-x^2} \, \exp[-\xi(1-x^2)]$. See H. Albrect [*et al.*]{} \[The ARGUS Collaboration\], Phys. Lett. [**B241**]{}, 278 (1990). K. Nakamura [*et al.*]{} \[The Particle Data Group\], J. Phys. [**G37**]{}, 075021 (2010). Z. Bai [*et al.*]{} \[The Mark III Collaboration\], Phys. Rev. Lett, [**65**]{}, 1309 (1990). M. Ablikim [*et al.*]{} \[The BES Collaboration\], Phys. Lett. [**B662**]{}, 330 (2008). M. Pivk and F. R. Le Diberder, Nucl. Instrum. Meth. [**A555**]{}, 356 (2005).
