Dalitz plot analysis of the decay B+- ->K+-K+-K-+

We analyze the three-body charmless decay B+- ->K+-K+-K-+ using a sample of 226.0 +- 2.5 million BBbar pairs collected by the BABAR detector. We measure the total branching fraction and CP asymmetry to be B = (35.2 +- 0.9 +-1.6) x 10^{-6} and A_CP = (-1.7 +- 2.6 +- 1.5)%. We fit the Dalitz plot distribution using an isobar model and measure the magnitudes and phases of the decay coefficients. We find no evidence of CP violation for the individual components of the isobar model. The decay dynamics is dominated by the K+K- S-wave, for which we perform a partial-wave analysis in the region m(K+K-)<2 GeV/c^2. Significant production of the f0(980) resonance, and of a spin zero state near 1.55 GeV/c^2 are required in the isobar model description of the data. The partial-wave analysis supports this observation.

We analyze the three-body charmless decay B ± → K ± K ± K ∓ using a sample of 226.0 ± 2.5 million BB pairs collected by the BABAR detector. We measure the total branching fraction and CP asymmetry to be B = (35.2 ± 0.9 ± 1.6) × 10 −6 and ACP = (−1.7 ± 2.6 ± 1.5)%. We fit the Dalitz plot distribution using an isobar model and measure the magnitudes and phases of the decay coefficients. We find no evidence of CP violation for the individual components of the isobar model. The decay dynamics is dominated by the K + K − S-wave, for which we perform a partial-wave analysis in the region m(K + K − ) < 2 GeV/c 2 . Significant production of the f0(980) resonance, and of a spin zero state near 1.55 GeV/c 2 are required in the isobar model description of the data. The partial-wave analysis supports this observation.

I. INTRODUCTION
Charmless decays of B mesons provide a rich laboratory for studying different aspects of weak and strong interactions. With recent theoretical progress in understanding the strong interaction effects, specific predictions for two-body pseudoscalar-pseudoscalar, B → P P, and pseudoscalar-vector, B → P V , branching fractions and asymmetries are available [1,2,3,4] and global fits to experimental data have been performed [5,6]. Improved experimental measurements of a comprehensive set of charmless B decays coupled with further theoretical progress hold the potential to provide significant constraints on the CKM matrix parameters and to discover hints of physics beyond the Standard Model in penguinmediated b → s transitions.
We analyze the decay B ± → K ± K ± K ∓ , dominated by the b → s penguin-loop transition, using 226.0 ± 2.5 million BB pairs collected by the BABAR detector [7] at the SLAC PEP-II asymmetric-energy B factory [8] operating at the Υ (4S) resonance. BABAR has previously measured the total branching fraction and asymmetry in this mode [9] and the two-body branching fractions B ± → K ± φ(1020), B ± → K ± χ c0 [10,11]. A comprehensive Dalitz plot analysis of B ± → K ± K ± K ∓ has been published by the Belle collaboration [12].

II. EVENT SELECTION
We consider events with at least four reliably reconstructed charged-particle tracks consistent with having originated from the interaction point. All three tracks forming a B ± → K ± K ± K ∓ decay candidate are required to be consistent with a kaon hypothesis using a particle identification algorithm that has an average efficiency of 94% within the acceptance of the detector and an average pion-to-kaon misidentification probability of 6%. * Also at Laboratoire de Physique Corpusculaire, Clermont-Ferrand, France † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università della Basilicata, Potenza, Italy We use two kinematic variables to identify the signal. The first is ∆E = E − √ s 0 /2, the difference between the reconstructed B candidate energy and half the energy of the e + e − initial state, both in the e + e − center-of-mass (CM) frame. For signal events the ∆E distribution peaks near zero with a resolution of 21 MeV. We require the candidates to have |∆E| < 40 MeV. The second variable is the energy-substituted mass where p B is the momentum of the B candidate and (E 0 , p 0 ) is the fourmomentum of the e + e − initial state, both in the laboratory frame. For signal events the m ES distribution peaks near the B mass with a resolution of 2.6 MeV/c 2 . We define a signal region (SR) with m ES ∈ (5.27, 5.29) GeV/c 2 and a sideband (SB) with m ES ∈ (5.20, 5.25) GeV/c 2 .
The dominant background is due to events from lightquark or charm continuum production, e + e − → qq, whose jet-like event topology is different from the more spherical B decays. We suppress this continuum background by requiring the absolute value of the cosine of the angle between the thrust axes of the B candidate and the rest of the event in the CM frame to be smaller than 0.95. Further suppression is achieved using a neural network with four inputs computed in the CM frame: the cosine of the angle between the direction of the B candidate and the beam direction; the absolute value of the cosine of the angle between the candidate thrust axis and the beam direction; and momentum-weighted sums over tracks and neutral clusters not belonging to the candidate, i p i and i p i cos 2 θ i , where p i is the track momentum and θ i is the angle between the track momentum direction and the candidate thrust axis. Figure 1 shows the m ES distribution of the 9870 events thus selected. The histogram is fitted with a sum of a Gaussian distribution and a background function having a probability density, √ s 0 and ξ is a shape parameter [13]. The binned maximum likelihood fit gives χ 2 = 104 for 100 bins and ξ = 21.1 ± 1.6. The ratio of the integrals of the background function over the signal region and the sideband yields an extrapolation coefficient R qq = 0.231 ± 0.007. The expected number of qq background events in the signal region is n SR qq = R qq (n SB − n SB BB ) = 972 ± 34, where n SB = 4659 is the number of events in the sideband from which we subtract the number of nonsignal BB background events n SB BB = 431 ± 19, estimated using a large number of simulated exclusive B decays. The expected number of BB background events in the signal region is n SR BB = 276±20, with B ± → DK ± decays giving the largest contribution.
We use a kinematic fit, constraining the mass of the selected candidates to the mass of the B meson. The threebody decay kinematics is described by two independent di-kaon invariant mass variables (m 2 23 , m 2 13 ) = (s 23 , s 13 ), where we order the same-sign kaons such that s 23 ≤ s 13 . The signal region contains 1769 B + and 1730 B − candidates whose Dalitz plot distribution is shown in Fig. 2.

III. ISOBAR MODEL FIT
We perform an extended binned maximum likelihood fit to the event distribution in Fig. 2 by binning the folded Dalitz plot into 292 non-uniform rectangular bins and minimizing the log of the Poisson likelihood ratio, is the number of observed signal region events in the i-th bin, assumed to be sampled from a Poisson distribution with mean µ i . In the limit of large statistics, the χ 2 LLR function has a χ 2 distribution and can be used to evaluate the goodness-of-fit. The expected number of events in the i-th bin is modeled as where the first term is the expected signal contribution, given by twice the bin integral of the square of the matrix element multiplied by the signal selection efficiency ǫ, determined as a function of (m 23 , m 13 ) using simulated signal events. The integral is multiplied by two because we use a folded Dalitz plot. We use the isobar model formalism [14,15] and describe the matrix element M as a sum of coherent contributions, M = N k=1 M k . The individual contributions are symmetrized with respect to the interchange of same-sign kaons, 1 ↔ 2, and are given by where ρ k e iφ k is a complex-valued decay coefficient, A k is the amplitude describing a K + K − system in a state with angular momentum J and invariant mass √ s 23 , P J is the Legendre polynomial of order J, and the helicity angle θ 13 between the direction of the bachelor recoil kaon 1 and kaon 3 is measured in the rest frame of kaons 2 and 3.
The model includes contributions from the φ(1020) and χ c0 intermediate resonances, which are clearly visible in Fig. 2. Following Ref. [12], we introduce a broad scalar resonance, whose interference with a slowly varying nonresonant component is used to describe the rapid decrease in event density around m(K + K − ) = 1.6 GeV/c 2 . Evidence of a possible resonant S-wave contribution in this region has been reported previously [16,17], however its attribution is uncertain: the f 0 (1370) and f 0 (1500) resonances are known to couple more strongly to ππ than to KK [18]; possible interpretations in terms of those states [19] must account for the fact that no strong [12,20]. The contribution of the f 0 (1710) resonance is included in the fit as a separate component and is found to be small. In the following, we designate the broad scalar resonance X 0 (1550) and determine its mass and width directly from the fit.
The contribution from a spin J resonance with mass m 0 and total width Γ 0 is described by a relativistic Breit-Wigner amplitude: . (

3)
F J is the Blatt-Weisskopf centrifugal barrier factor [21] for angular momentum J: h , and R represents the effective radius of the interaction volume for the resonance; we use R = 4.0 GeV −1 (0.8 fm) [22]. In the formulation of Eq. (3), only the centrifugal barrier factor for the decay of a spin J resonance into two pseudoscalar kaons is included; we have ignored the corresponding centrifugal barrier factor for the two-body decay of a B meson into a pseudoscalar kaon and a spin J resonance. The effect of this approximation on the parameterization of B ± → K ± φ(1020), the only component with J > 0, is negligible. Unless otherwise specified, all resonance parameters are taken from Ref. [18]. The term ∆Γ(s), parameterizing the mass dependence of the total width, is in general given by ∆Γ(s) = i ∆Γ i (s), where the sum is over all decay modes of the resonance, and ∆Γ i (m 2 0 ) ≡ 0. The χ c0 has many decay modes, the decay modes of the f 0 (1710) are not well established, and decay modes other than K + K − of the possible X 0 (1550) resonance are unknown; in all these cases we set ∆Γ(s) = 0 and neglect the mass dependence of the total width. For the φ(1020) resonance we use ∆Γ(s) = ∆Γ 1 (s) + ∆Γ 2 (s), where , and the mass dependence of the partial width for the two-body vector to pseudoscalar-pseudoscalar decay φ(1020) → hh is parameterized as [12,20], and a recent measurement of g K /g π , the ratio of the f 0 (980) coupling constants to KK and ππ [23], motivate us to include an f 0 (980) contribution using a coupled-channel amplitude parameterization: where ̺ π = 2/3 1 − 4m 2 π ± /s + 1/3 1 − 4m 2 π 0 /s, ̺ K = 1/2 1 − 4m 2 K ± /s + 1/2 1 − 4m 2 K 0 /s, and we use g K /g π = 4.21 ± 0.25 ± 0.21, m 0 = 0.965 ± 0.008 ± 0.006 GeV/c 2 and g π = 0.165 ± 0.010 ± 0.015 GeV/c 2 [23].
We have investigated two theoretical models of the nonresonant component [24,25] and found that neither describes the data adequately. We therefore include an Swave nonresonant component expanded beyond the usual constant term as In the following we fix β = 0, and incorporate the M NR contribution over the entire Dalitz plot, thus effectively employing the same parameterization as in Ref. [12]. We fit for the magnitudes and phases of the decay coefficients, the mass and width of the X 0 (1550), and the nonresonant component shape parameter α. As the overall complex phase of the isobar model amplitude is arbi-trary, we fix the phase of the nonresonant contribution to zero, leaving 14 free parameters in the fit. The number of degrees of freedom is 292 − 14 = 278. We perform multiple minimizations with different starting points and find multiple solutions clustered in pairs, where the solutions within each pair are very similar, except for the magnitude and phase of the χ c0 decay coefficient. The twofold ambiguity arises from the interference between the narrow χ c0 and the nonresonant component, which is approximately constant across the resonance. The highest-likelihood pair has χ 2 LLR = (346.4, 352.0); the second best pair has χ 2 LLR = (362.4, 368.7). The least significant components are the f 0 (980) and the f 0 (1710). Their omission from the fit model degrades the best fit from χ 2 LLR = 346.4 to 363.9 and 360.7, respectively. The invariant-mass projections of the best fit are shown in Fig. 3. The goodness-of-fit is χ 2 = 56 for 56 bins in the m 23 projection and χ 2 = 66 for 63 bins in the m 13 projection. The sharp peak in the BB background distribution in the m 23 projection is due to the contribution from the B ± → DK ± backgrounds. The fit gives α = 0.152 ± 0.011 GeV −2 c 4 , m 0 (X 0 ) = 1.539 ± 0.020 GeV/c 2 , and Γ 0 (X 0 ) = 0.257 ± 0.033 GeV/c 2 . The fitted values of the shape parameter α and the resonance mass are consistent with the values in Ref. [12], but our preferred value for the width is significantly larger. The results of the best isobar model fit are summarized in Table I, where we have also included the results for the χ c0 component from the second solution in the highest-likelihood pair, labeled χ II c0 , and component fit fractions, where the integrals are taken over the entire Dalitz plot. The sum of the component fit fractions is significantly larger than one due to large negative interference in the scalar sector [26].

IV. BRANCHING FRACTIONS AND ASYMMETRIES
To search for possible direct CP violation we extend the isobar model formalism by defining charge-dependent decay coefficients: where A k is the CP asymmetry of the k-th component, and δφ k = φ − k − φ + k . We modify the likelihood to be the product of the likelihoods for the two charges and repeat the fit. The phase of the nonresonant component is fixed to zero for both charges. The results are given in the last three columns of Table I in terms of the fitted CP asymmetry values, the symmetric 90% confidence level intervals around them, and the fitted phase differences between the charge-dependent decay coefficients.
The asymmetry intervals are estimated by fitting Monte Carlo simulated samples generated according to the parameterized model of the nominal asymmetry fit. There is no evidence of statistically significant CP violation for any of the components.
Taking into account the signal Dalitz plot distribution, as described by the isobar model fit, the average signal efficiency isε = 0.282±0.011, where the uncertainty is evaluated using control data samples, and is primarily due to the uncertainties in tracking and particle identification efficiencies. The total branching fraction is B(B ± → K ± K ± K ∓ ) = (35.2 ± 0.9 ± 1.6) × 10 −6 , where the first error is statistical and the second is systematic. The fit fraction of the isobar model terms that do not involve the χ c0 resonance is (95.0±0.6±1.1)% for the best fit, giving B(B ± → K ± K ± K ∓ ) = (33.5 ± 0.9 ± 1.6)× 10 −6 if intrinsic charm contributions are excluded. The total asymmetry is The systematic error for the overall branching fraction is obtained by combining in quadrature the 3.9% efficiency uncertainty, a 1.1% uncertainty on the total number of B + B − pairs, a 0.7% uncertainty due to the modeling of BB backgrounds, and a 1.4% uncertainty due to the uncertainty arising from the uncertainty on the R qq sideband extrapolation coefficient. The 1.5% systematic uncertainty for the asymmetry is due to possible charge asymmetry in kaon tracking and particle identification efficiencies, evaluated using data control samples. Where appropriate, the systematic uncertainties discussed above have been propagated to estimate the uncertainties on the leading isobar model fit results. We have also evaluated the systematic uncertainties due to the parameterization of resonance lineshapes by varying the parameters of all resonances within their respective uncertainties. Uncertainties arising from the distortion of narrow resonance lineshapes due to finite detector resolution and, for candidates containing a φ(1020) resonance produced in the qq continuum, due to the kinematic fit, have also been studied.
The partial branching fractions for B ± → K ± f 0 (980) measured in the K ± K ± K ∓ and K ± π ± π ∓ final states are related by the ratio where 3/4 is an isospin factor, and I K /I π is the ratio of the integrals of the square of the f 0 (980) amplitude The magnitudes and phases of the decay coefficients, fit fractions, two-body branching fractions, CP asymmetries, symmetric 90% confidence level CP asymmetry intervals around the nominal value, and the phase differences between the charge-dependent decay coefficients for the individual components of the isobar model fit.

V. PARTIAL-WAVE ANALYSES
We further study the nature of the dominant S-wave component by considering the interference between the low-mass and the high-mass scattering amplitudes in the region m 23 ∈ (1.1, 1.8) GeV/c 2 , m 13 > 2 GeV/c 2 . The matrix element is modeled as where ρ S and φ S refer to the S-wave and are taken to be constant within each bin of the s 23 variable and the nonresonant amplitude parameterization is taken from the fit to the high-mass region. The partial-wave expansion truncated at the S-wave describes the data adequately; the magnitude of the S-wave in each bin is readily determined. Because of the mass dependence of the nonresonant component, the phase of the S-wave can also be determined, albeit with a sign ambiguity and rather large errors for bins with a small number of entries or small net variation of the nonresonant component. The results are shown in Fig. 4, with the S-wave component of the isobar model fit overlaid for comparison. Continuity requirements allow us to identify two possible solutions for the phase; the solution labeled by black squares is consistent with a rapid counterclockwise motion in the Argand plot around m(K + K − ) = 1.55 GeV/c 2 , which is accommodated in the isobar model as the contribution of the X 0 (1550).
Isospin symmetry relates the measurements in B ± → K ± K ± K ∓ and B 0 → K + K − K 0 S [29]. Our results for the K + K − S-wave can therefore be used to estimate a potentially significant source of uncertainty in the measurements of sin2β in B 0 → φ(1020)K 0 S [30,31] due to the   contribution of a CP -even S-wave amplitude. We perform a partial-wave analysis in the region m 23 (K + K − ) ∈ (1.013, 1.027) GeV/c 2 , which we assume to be dominated by the low-mass P -wave, due to the contribution of the φ(1020) resonance, and a low-mass S-wave. The matrix element is modeled as where the low-mass S-wave is taken to be constant over the small s 23 interval considered. The fit results for the P -wave are shown in Fig. 5, with a Breit-Wigner fit of the φ(1020) resonance overlaid for comparison. For the Swave we get ρ 2 S = (3.4 ± 2.5) × 10 2 GeV −4 c 8 and compute its fraction in this region using Eq. (7) to be (9 ± 6)%.
We also consider the region 2m K + < m(K + K − ) < 1.006 GeV/c 2 , in the immediate vicinity of the K + K − threshold. The contribution of the φ(1020) resonance tail in this region is suppressed by the centrifugal barrier and is estimated to be smaller than 10%. We fit ρ 2 S = (6.1 ± 1.6) × 10 2 GeV −4 c 8 for the magnitude of the Swave in this region. The fits in the vicinity of the K + K − threshold and in the region around the φ(1020) resonance indicate a threshold enhancement of the S-wave, which is accommodated in the isobar model by the contribution of the f 0 (980) resonance as shown in the inset of Fig. 4.

VI. CONCLUSIONS
In conclusion, we have measured the total branching fraction and the CP asymmetry in B ± → K ± K ± K ∓ . An isobar model Dalitz plot fit and a partial-wave analysis of the K + K − S-wave show evidence of large contributions from a broad X 0 (1550) scalar resonance, a massdependent nonresonant component, and an f 0 (980) resonance. The ratio of B ± → K ± f 0 (980) two-body branching fractions measured by BABAR in B ± → K ± π ± π ∓ and B ± → K ± K ± K ∓ is consistent with the measurement of g K /g π by the BES collaboration, albeit with large errors. Our isobar model fit results are substantially different from those obtained in Ref. [12] due to the larger fitted width of the X 0 (1550) and the inclusion of the f 0 (980) component in the isobar model. Our results for the B(B ± → K ± φ(1020)) and B(B ± → K ± χ c0 ) branching fractions are in agreement with the previous results from BABAR [10,11], which they supersede, and from other experiments [12,27,28]. We have measured the CP asymmetries and the phase differences between the charge-dependent decay coefficients for the individual components of the isobar model and found no evidence of direct CP violation.
We wish to dedicate this paper to Prof. Richard E. (Dick) Dalitz, inventor of the "Phase Space Plot" (as he called it) with which his name is so intimately linked, in tribute as much to his qualities as a gentleman as to his gifts as a physicist.