Search for rare B(s)0→μ+μ -μ+μ- decays

A search for the decays Bs ! þ þ and B ! þ þ is performed using data, corresponding to an integrated luminosity of 1:0 fb , collected with the LHCb detector in 2011. The number of candidates observed is consistent with the expected background and, assuming phase-space models of the decays, limits on the branching fractions are set: BðBs! þ þ Þ<1:6ð1:2Þ 10 8 and BðB0 ! þ þ Þ< 6:6ð5:3Þ 10 9 at 95% (90%) confidence level. In addition, limits are set in the context of a supersymmetric model which allows for the B0ðsÞ meson to decay into a scalar (S) and pseudoscalar particle (P), where S and P have masses of 2:5 GeV=c and 214:3 MeV=c, respectively, both resonances decay into þ . The branching fraction limits for these decays are BðBs ! SPÞ< 1:6ð1:2Þ 10 8 and BðB0 ! SPÞ< 6:3ð5:1Þ 10 9 at 95% (90%) confidence level.

R. Aaij et al. * (LHCb Collaboration) (Received 6 March 2013;published 21 May 2013) A search for the decays B 0 s ! þ À þ À and B 0 ! þ À þ À is performed using data, corresponding to an integrated luminosity of 1:0 fb À1 , collected with the LHCb detector in 2011. The number of candidates observed is consistent with the expected background and, assuming phase-space models of the decays, limits on the branching fractions are set: BðB 0 s ! þ À þ À Þ<1:6ð1:2ÞÂ10 À8 and BðB 0 ! þ À þ À Þ < 6:6ð5:3Þ Â 10 À9 at 95% (90%) confidence level. In addition, limits are set in the context of a supersymmetric model which allows for the B 0 ðsÞ meson to decay into a scalar (S) and pseudoscalar particle (P), where S and P have masses of 2:5 GeV=c and 214:3 MeV=c, respectively, both resonances decay into þ À . The branching fraction limits for these decays are BðB 0 s ! SPÞ < 1:6ð1:2Þ Â 10 À8 and BðB 0 ! SPÞ < 6:3ð5:1Þ Â 10 À9 at 95% (90%) confidence level. DOI: 10.1103/PhysRevLett.110.211801 PACS numbers: 13.20.He, 12.15.Mm, 12.60.Jv The decays B 0 ðsÞ ! þ À þ À are strongly suppressed in the standard model (SM). Significant enhancements in the branching fractions can occur in beyond the SM theories [1,2]. For example, in minimal supersymmetric models (MSSM) the decay is mediated via the new scalar S, and pseudoscalar P sgoldstino particles, which both decay into þ À . The sgoldstinos can couple to SM particles via type-I couplings, where one sgoldstino couples to two SM fermions and via type-II couplings, where both S and P couple to two SM fermions in a four-prong vertex [1]. Such searches are also of particular interest since the HyperCP Collaboration found evidence of the decay AE þ ! p þ À , which is consistent with the decay AE þ ! pP and P ! þ À with a P particle mass of 214:3 AE 0:5 MeV=c 2 [3]. The inclusion of charge conjugated processes is implied throughout this Letter. The Belle Collaboration has reported a search for a P particle in B 0 ! VPð! þ À Þ decays, where V is either a K Ã0 ð892Þ or a 0 ð770Þ meson [4]. The corresponding upper limits constrain the type-I couplings of sgoldstinos to the SM particles. The search presented here is sensitive to the decays B 0 ðsÞ ! SP in which both the S and P particle decay vertices are not significantly displaced from the B 0 ðsÞ decay vertex and the P particle has a similar mass to that reported by the HyperCP collaboration. Such final states probe type-II couplings, which might dominate for general low energy supersymmetry breaking models and have not been probed by previous experiments. The present search is therefore complementary to the Belle study. Moreover, the four-muon final state is essentially background free and therefore ideally suited to search for new physics signatures. The present search is also sensitive to B 0 ðsÞ ! þ À þ À decays which are not propagated through intermediate resonant structures.
The dominant SM decay of a B meson into a four-muon final state is B 0 s ! J=c ð1020Þ, where both the J=c and the mesons decay into two muons [ Fig. 1(a)]. In this Letter, this is referred to as the resonant decay mode. The branching fraction for B 0 s !J=c ð! þ À Þð! þ À Þ is calculated as the product of BðB 0 s ! J=c Þ, BðJ=c ! þ À Þ and Bð ! þ À Þ [5], resulting in a value of ð2:3 AE 0:9Þ Â 10 À8 . The main SM nonresonant B 0 ðsÞ ! þ À þ À decay mode is B 0 ðsÞ ! þ À ð! þ À Þ, where one muon pair is produced via an electroweak penguin or box diagram and the other via a virtual photon [ Fig. 1(b)]. The branching fraction of B 0 ðsÞ ! þ À ð! þ À Þ is expected to be less than 10 À10 [6]. The diagram for the MSSM decay mode B 0 ðsÞ !SP is shown in Fig. 1(c).
This Letter presents a search for the decays B 0 ðsÞ ! þ À þ À using data corresponding to an integrated luminosity of 1:0 fb À1 of pp collision data at ffiffi ffi s p ¼ 7 TeV collected with the LHCb detector in 2011. The resonant B 0 s ! J=c decay mode is removed in the signal selection and is used as a control channel to develop the selection criteria. The decay B 0 ! J=c K Ã0 , where J=c ! þ À and K Ã0 ! K þ À , is used as a normalization channel to measure the branching fractions of the B 0 ðsÞ ! þ À þ À decays. The detector [7] is a single-arm forward spectrometer, covering the pseudorapidity range 2 < < 5, designed for the study of particles containing b or c quarks. The detector includes a high precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet, and three stations of silicon-strip detectors and straw drift tubes placed downstream. The *Full author list given at the end of the article.
Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. magnet has a bending power of about 4 Tm. The combined tracking system has a momentum resolution Áp=p that varies from 0.4% at momenta of 5 GeV=c to 0.6% at 100 GeV=c. The impact parameter (IP) resolution is 20 m for tracks with high transverse momentum (p T ). Charged particles are identified using two ring-imaging Č erenkov detectors. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The trigger consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage that performs a full event reconstruction. Events are selected by the single-muon, dimuon, and generic b-hadron triggers described in Ref. [8].
Simulated events are generated using PYTHIA 6.4 [9] configured with the parameters detailed in Ref. [10]. Final state QED radiative corrections are included using the PHOTOS package [11]. The EVTGEN [12] and GEANT4 [13,14] packages are used to generate hadron decays and simulate interactions in the detector, respectively.
Signal B 0 ðsÞ ! þ À þ À candidates are selected by applying cuts on the final state muons and the reconstructed B 0 ðsÞ meson. Each final state muon candidate track is required to be matched with hits in the muon system [15]. The muon candidates are required to have particle identification (PID) criteria consistent with those of a muon and not those of a kaon or pion. This is determined by calculating an overall event likelihood for the distribution of Č erenkov photons detected by the ring-imaging Č erenkov system being consistent with a given particle hypothesis. To assess a particle hypotheses the difference between the logarithm of its likelihood and the pion hypothesis likelihood (DLL) is computed. Each muon candidate is required to have DLLðK ÀÞ<0 and DLLð À Þ > 0. The PID selection criteria yield a muon selection efficiency of 78.5%, the corresponding efficiency for mis-identifying a pion (kaon) as a muon is 1.4% (< 0:1%). All muon candidates are required to have a track fit 2 per degree of freedom of less than 5. Selection criteria are applied on the consistency of the muons to originate from a secondary vertex rather than a primary vertex. The muons are each required to have the difference between the 2 of the primary vertex formed with and without the considered tracks, 2 IP , to be greater than 16. The B 0 ðsÞ candidates are formed from two pairs of oppositely charged muons. The B 0 ðsÞ decay vertex 2 is required to be less than 30 to ensure that the four muons originate from a single vertex. In addition, the reconstructed B 0 ðsÞ meson is required to have a 2 IP less than 9 and hence be consistent with originating from a primary vertex.
The B 0 ðsÞ candidates are divided into two samples according to the invariant mass of the muon pairs. Signal nonresonant candidates are required to have all four þ À invariant mass combinations outside the respective and J=c mass windows of 950-1090 MeV=c 2 and 3000-3200 MeV=c 2 . The four-muon mass resolution is estimated by the simulation to be around 19 MeV=c 2 for both B 0 and B 0 s decay modes. The signal candidates are selected in a four-muon invariant mass window of AE40 MeV=c 2 around the world average B 0 ðsÞ mass [5]. Candidate B 0 s ! J=c decays are used to optimize the selection criteria described above. The B 0 s ! J=c candidates are selected by requiring the invariant mass of one muon pair to be within the mass window and that of the other pair to be within the J=c mass window.
The selection criteria are chosen by applying initial cut values that select generic B meson decays. These values are then further optimized using B 0 s ! J=c candidates. After applying the selection, seven B 0 s ! J=c candidates are observed in the signal sample. This is consistent with the expected B 0 s ! J=c yield, 5:5 AE 2:3, calculated by normalizing to the B 0 ! J=c K Ã0 decay mode.
The dominant B 0 ðsÞ ! þ À þ À background is combinatorial, where a candidate B 0 ðsÞ vertex is constructed from four muons that did not originate from a single B 0 ðsÞ meson. Sources of peaking background are estimated to be FIG. 1. Feynman diagrams for the B 0 s ! þ À þ À and B 0 ! þ À þ À decays, via (a) the resonant B 0 s ! J=c SM channel, (b) the nonresonant SM channel, and (c) the supersymmetric channel. The latter is propagated by scalar S and pseudoscalar P sgoldstino sfermions. negligible, the largest of these is due to B 0 ! c ð2SÞ Â ð! þ À ÞK Ã0 ð! K þ À Þ decays, which has an expected yield of 0:44 AE 0:06 events across the entire four-muon invariant mass range of 4776-5426 MeV=c 2 .
To evaluate the combinatorial background, a single exponential probability distribution function (PDF) is used to fit the events in mass ranges of 4776-5220 and 5426-5966 MeV=c 2 , where no signal is expected. Extrapolating the PDF into the B 0 ðB 0 s Þ signal window results in an expected background of 0:38 þ0:23 À0:17 (0:30 þ0:22 À0:20 ) events. Linear and double exponential fit models give consistent background yields.
The branching fraction of the B 0 ðsÞ ! þ À þ À decay is measured relative to that of the normalization channel B 0 ! J=c K Ã0 , this avoids uncertainties associated with the B 0 production cross section and the integrated luminosity. This normalization channel has the same topology as the signal channel and two muons as final state particles. The S-wave component from the nonresonant decay B 0 !J=c K þ À is removed from the present search. The branching fraction of B 0 !J=c ð! þ À ÞK Ã0 ð!K þ À Þ is calculated as the product of BðB 0 ! J=c K Ã0 Þ, with the S-wave component removed [16], BðK Ã0 ! K þ À Þ and BðJ=c ! þ À Þ [5]; the resulting branching fraction is ð5:10 AE 0:52Þ Â 10 À5 .
The muon PID and kinematic selection criteria for the normalization channel are identical to those applied to the signal channel. In addition to these criteria, the kaon is required to have its PID consistent with that of a kaon and not a pion, and vice versa for the pion. This removes B 0 ! J=c K Ã0 candidates where the kaon and pion mass hypotheses are exchanged. The K Ã0 ðJ=c Þ meson is selected by applying a mass window of AE100 MeV=c 2 ðAE50 MeV=c 2 Þ around the world average invariant mass [5].
There are four main sources of peaking backgrounds for B 0 ! J=c K Ã0 decays. The first is the decay B þ !J=c K þ combined with a pion from elsewhere in the event, which is removed by applying a veto on candidates with a K þ þ À invariant mass within AE60 MeV=c 2 of the world average B þ mass [5]. The second arises from B 0 s ! J=c decays, where the meson decays to two kaons, one of which is misidentified as a pion. This background is suppressed with a veto on candidates with a K þ À invariant mass, where the pion is assigned a kaon mass hypothesis, within AE70 MeV=c 2 of the world average mass [5]. The third source of background is from the decay " B 0 s ! J=c K Ã0 , which is included in the K þ À þ À invariant mass fit described below. The last source of background is the S-wave component of the decay mode B 0 ! J=c K þ À ; this is discussed later.
The yield of the normalization channel is determined by fitting the K þ À þ À invariant mass distribution with a combination of three PDFs, the B 0 ! J=c K Ã0 and the " B 0 s ! J=c K Ã0 signal PDFs consist of the sum of a Gaussian and a Crystal Ball function [17], centered around the respective world average B 0 and B 0 s masses and the background PDF consists of a single exponential. The resulting fit is shown in Fig. 2. The B 0 ! J=c K Ã0 mass peak has a Gaussian standard deviation width of 15:9 AE 0:6 MeV=c 2 and contains 31800 AE 200 candidates.
The branching fraction of the B 0 ðsÞ ! þ À þ À decay is calculated using ; ( where BðB 0 ! J=c K Ã0 Þ is the branching fraction of the normalization channel [5,16] and B 0 !J=c K Ã0 and B 0 ðsÞ ! þ À þ À are the efficiencies for triggering, reconstructing, and selecting the normalization and signal channel events, respectively. The efficiencies are calculated using simulated events and are cross-checked on data. The yields of the normalization and signal channels are N B 0 !J=c K Ã0 and N B 0 ðsÞ ! þ À þ À , respectively. The relative production fraction for B 0 and B 0 ðsÞ mesons, f dðsÞ =f d , is measured to be f s =f d ¼ 0:256 AE 0:020 for B 0 s decays [18] and taken as unity for B 0 decays. The factor accounts for the efficiency-corrected S-wave contribution to the normalisation channel yield; is calculated to be 1:09 AE 0:09, using the technique described in Ref. [19]. The PID components of the selection efficiencies are determined from data calibration samples of kaons, pions, and muons. The kaon and pion samples are obtained from D 0 ! K À þ decays, where the D 0 meson is produced via D Ãþ ! D 0 þ decays. The muon sample is obtained from B þ !J=c ð! þ À ÞK þ decays. The calibration samples are divided into bins of momentum, pseudorapidity, and the number of charged tracks in the event. This procedure corrects for differences between the kinematic and track multiplicity distributions of the simulated and the calibration event samples. Two models are used to simulate B 0 ðsÞ ! þ À þ À decays: (i) the phase space model, where the B 0 ðsÞ mass is fixed and the kinematics of the final state muons are distributed according to the available phase-space and (ii) the MSSM model, which describes the decay mode B 0 ðsÞ ! SP. In the MSSM model the pseudoscalar particle P is a sgoldstino of mass 214:3 MeV=c 2 , consistent with results from the HyperCP experiment [3]. The decay widths of S and P are set to 0:1 MeV=c 2 . The scalar sgoldstino S mass is set to 2:5 GeV=c 2 . If the mass of S is varied across the allowed phase space of the B 0 ðsÞ ! SP decay, the relative change in B 0 ðsÞ !SP from the central values varies from À12:6% to þ17:2%.
The calculated efficiencies of all the simulated decay modes are shown in Table I. The total efficiencies of the MSSM models are comparable to those for the phase space models, indicating that the present search has approximately the same sensitivity to new physics models, which feature low mass resonances, as to the phase space models.
Systematic uncertainties enter into the calculation of the limits on BðB 0 ðsÞ ! þ À þ À Þ through the various elements of Eq. (1). The largest uncertainty arises from the branching fraction of B 0 ! J=c K Ã0 , which is known to a precision of 10.2% [5,16]. An uncertainty arises due to the correction for the B 0 ! J=c K þ À S-wave contribution. This is conservatively estimated to be 8.3%, which is the maximum relative change in when it is calculated: (i) by using the angular acceptance from simulated events and (ii) by performing a fit with the physics parameters of the decay fixed and the angular acceptance parameterised, the coefficients of which are left free in the fit. An uncertainty of 7.8% is introduced in the calculation of BðB 0 s ! þ À þ À Þ, due to the uncertainty on f s =f d [18]. The 4.4% systematic uncertainty associated with the trigger efficiency is calculated as the relative difference between the data and simulation efficiencies of the trigger selection criteria applied to the normalization channel. The efficiency in data is calculated using the method described in Ref. [20]. Small differences are seen between the data and the simulated events for the track 2 IP distributions and the efficiency for reconstructing individual tracks. The distributions of these quantities are corrected in the simulation to resemble the data using data-driven methods and the associated uncertainty is assessed by varying the magnitude and the configuration of the corrections. The relative systematic uncertainty assigned to the ratio of efficiencies, B 0 !J=c K Ã0 = B 0 s ! þ À þ À , is calculated to be 5.2%. The 4.1% uncertainty on the PID selection efficiency is the maximum relative change in B 0 !J=c K Ã0 = B 0 s ! þ À þ À that results from applying different binning schemes to the PID calibration samples. This uncertainty includes effects associated with muon, kaon and pion identification. The statistical uncertainty associated with the size of the simulated event samples is 1.3% for both B 0 and B 0 s modes. The uncertainty on the B 0 ! J=c K Ã0 yield is 0.6%. Table II summarizes the systematic uncertainties for both the BðB 0 s ! þ À þ À Þ and BðB 0 ! þ À þ À Þ branching fractions. The combined systematic uncertainty for the B 0 s and B 0 modes is 17.2% and 15.4%, respectively. The same uncertainties apply for BðB 0 s ! SPÞ and BðB 0 ! SPÞ. The nonresonant four-muon invariant mass range and the background fit are shown in Fig. 3. One event is observed in the B 0 signal window and zero events are observed in the B 0 s window. These observations are consistent with the expected background yields. The CL s method [21,22] is used to set upper limits on the branching fractions. The 95% (90%) confidence level limits for the nonresonant B 0 ðsÞ ! þ À þ À decay modes in the phase space model are BðB 0 s ! þ À þ À Þ < 1:6ð1:2Þ Â 10 À8 ; BðB 0 ! þ À þ À Þ < 6:6ð5:3Þ Â 10 À9 : The corresponding limits for the MSSM model with B 0 ðsÞ ! SP and the mass of PðSÞ set to 214:3 MeV=c 2 ð2:5 GeV=c 2 Þ, are  BðB 0 s ! SPÞ < 1:6ð1:2Þ Â 10 À8 ; BðB 0 ! SPÞ < 6:3ð5:1Þ Â 10 À9 : Varying the mass of S across the allowed phase space of the B 0 ðsÞ ! SP decay, from 211 MeV=c 2 to 5065 ð5153Þ MeV=c 2 for B 0 ðB 0 s Þ, results in a relative change in the 95% confidence level limit from À23% to þ6% for both B 0 and B 0 s decay modes. In summary, a search for the decays B 0 ðsÞ ! þ À þ À has been performed and first limits on the branching fractions for these decay modes have been set. These limits probe the upper regions of the parameter space of the B 0 ðsÞ ! SP decay, and, in particular, set the first constraints on type-II couplings [1].
We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq, FAPERJ and FINEP (Brazil); NSFC (China); CNRS/IN2P3 and Region Auvergne (France); BMBF, DFG, HGF and MPG (Germany); SFI (Ireland); INFN (Italy); FOM and NWO (The Netherlands); SCSR (Poland); ANCS/IFA (Romania); MinES, Rosatom, RFBR and NRC ''Kurchatov Institute'' (Russia); MinECo, XuntaGal and GENCAT (Spain); SNSF and SER (Switzerland); NAS Ukraine (Ukraine); STFC (United Kingdom); NSF (USA). We also acknowledge the support received from the ERC under FP7. The Tier1 computing centres are supported by IN2P3 (France), KIT and BMBF (Germany), INFN (Italy), NWO and SURF (Netherlands), PIC (Spain), GridPP (United Kingdom). We are thankful for the computing resources put at our disposal by Yandex LLC (Russia), as well as to the communities behind the multiple open source software packages that we depend on.