Search for rare $B^0_{(s)}\rightarrow \mu^+ \mu^- \mu^+ \mu^-$ decays

\noindent A search for the decays $B^0_{s}\rightarrow \mu^+ \mu^- \mu^+ \mu^-$ and $B^0 \rightarrow \mu^+ \mu^- \mu^+ \mu^-$ is performed using data, corresponding to an integrated luminosity of 1.0\ensuremath{\{\,fb}^{-1}}\xspace, 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: \{${\ensuremath{\cal B}\xspace}(B^0_{s}\rightarrow \mu^+ \mu^- \mu^+ \mu^-)<1.6 \ (1.2) \times 10^{-8}$} and \{${\ensuremath{\cal B}\xspace}(B^0 \rightarrow \mu^+ \mu^- \mu^+ \mu^-)<6.6 \ (5.3) \times 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 and 214.3 MeV, respectively, both resonances decay into $\mu^+\mu^-$. The branching fraction limits for these decays are \{${\ensuremath{\cal B}\xspace}(\ensuremath{B^0_{s}\rightarrow SP}\xspace)<1.6 \ (1.2) \times 10^{-8}$} and \{${\ensuremath{\cal B}\xspace}(\ensuremath{B^0\rightarrow SP}\xspace)<6.3 \ (5.1) \times 10^{-9}$} at 95% (90%) confidence level.

a Also at P.N.Lebedev Physical Institute, Russian Academy of Science (LPI RAS), Moscow, Russia b Also at Università di Bari, Bari, Italy c Also at Università di Bologna, Bologna, Italy d Also at Università di Cagliari, Cagliari, Italy e Also at Università di Ferrara, Ferrara, Italy f Also at Università di Firenze, Firenze, Italy g Also at Università di Urbino, Urbino, Italy h Also at Università di Modena e Reggio Emilia, Modena, Italy i Also at Università di Genova, Genova, Italy j Also at Università di Milano Bicocca, Milano, Italy k Also at Università di Roma Tor Vergata, Roma, Italy l Also at Università di Roma La Sapienza, Roma, Italy m Also at Università della Basilicata, Potenza, Italy n Also at LIFAELS, La Salle, Universitat Ramon Llull, Barcelona, Spain o Also at Hanoi University of Science, Hanoi, Viet Nam p Also at Massachusetts Institute of Technology, Cambridge, MA, United States 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 Σ + → pµ + µ − , which is consistent with the decay Σ + → pP and P → µ + µ − with a P particle mass of 214.3 ± 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 → V P (→ µ + µ − ) 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/ψφ(1020), where both the J/ψ 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 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 in-tegrated luminosity of 1.0 fb −1 of pp collision data at √ s = 7 TeV collected with the LHCb detector in 2011.The resonant B 0 s → J/ψφ 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/ψK * 0 , where J/ψ → µ + µ − 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 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 re- the resonant B 0 s → J/ψφ SM channel, (b) the nonresonant SM channel and (c) the supersymmetric channel.The latter is propagated by scalar S and pseudoscalar P sgoldstino sfermions.quired 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 chi 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/ψ 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 ±40 MeV/c 2 around the world average B 0 (s) mass [5].Candidate B 0 s → J/ψφ decays are used to optimize the selection criteria described above.The B 0 s → J/ψφ 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/ψ 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/ψφ candidates.
After applying the selection, seven B 0 s → J/ψφ candidates are observed in the signal sample.This is consistent with the expected B 0 s → J/ψφ yield, 5.5 ± 2.3, calculated by normalizing to the B 0 → J/ψ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 negligible, the largest of these is due to B 0 → ψ(2S) (→ µ + µ − ) K * 0 (→ K + π − ) decays, which has an expected yield of 0.44 ± 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 muon PID and kinematic selection criteria for the normalisation 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/ψK * 0 candidates where the kaon and pion mass hypotheses are exchanged.The K * 0 (J/ψ) meson is selected by applying a mass window of ±100 MeV/c 2 (±50 MeV/c 2 ) around the world average invariant mass [5].
There are four main sources of peaking backgrounds for B 0 → J/ψK * 0 decays.The first is the decay B + → J/ψ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 ±60 MeV/c 2 of the world average B + mass [5].The second arises from B 0 s → J/ψφ 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 ±70 MeV/c 2 of the world average φ mass [5].The third source of background is from the decay B 0 s → J/ψ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/ψK + π − , this is discussed later.
The yield of the normalisation channel is determined by fitting the K + π − µ + µ − invariant mass distribution with a combination of three PDFs, the B 0 → J/ψK * 0 and the B 0 s → J/ψK * 0 signal PDFs consist of the sum of a Gaussian and a Crystal Ball function [17], centred 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/ψK * 0 mass peak has a Gaussian standard deviation width of 15.9 ± 0.6 MeV/c 2 and contains 31 800 ± 200 candidates.
The branching fraction of the B 0 (s) → µ + µ − µ + µ − decay is calculated using where B(B 0 → J/ψK * 0 ) is the branching fraction of the normalisation channel [5,16] and B 0 →J/ψK * 0 and B 0 (s) →µ + µ − µ + µ − are the efficiencies for triggering, reconstructing and selecting the normalisation and signal channel events, respectively.The efficiencies are calculated using simulated events and are cross-checked on data.The yields of the normalisation and signal channels are N B 0 →J/ψ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 ± 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 ± 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/ψ (→ µ + µ − )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 1.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/ψ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/ψ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 associ- (s) → µ + µ − µ + µ − candidates.The solid (dashed) black lines indicate the boundaries of the B 0 s (B 0 ) signal window.The blue curve shows the model used to fit the mass sidebands and extract the expected number of combinatorial background events in the B 0 s and B 0 signal regions.Only events in the region in which the line is solid have been considered in the fit.ated with the trigger efficiency is calculated as the relative difference between the data and simulation efficiencies of the trigger selection criteria applied to the normalisation 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 datadriven 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/ψ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/ψ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/ψK * 0 yield is 0.6 %.Table 2 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 ).

Figure 2 :
Figure 2: (color online).Invariant mass distribution of K + π − µ + µ − candidates.The B 0 and B 0 s signal distributions are shown by short-dashed black and long-dashed red (gray) lines, respectively.The background shape is denoted by a long-dashed black line.The total fit result is shown as a solid blue (black) line.The inset shows the mass distribution centered around the B 0 s mass.

Figure 3 :
Figure 3: Invariant mass distribution of nonresonant B 0(s) → µ + µ − µ + µ − candidates.The solid (dashed) black lines indicate the boundaries of the B 0 s (B 0 ) signal window.The blue curve shows the model used to fit the mass sidebands and extract the expected number of combinatorial background events in the B 0 s and B 0 signal regions.Only events in the region in which the line is solid have been considered in the fit.

Table 1 :
Combined reconstruction and selection efficiencies of all the decay modes considered in the analysis.The uncertainties shown are statistical.

Table 2 :
Systematic uncertainties on the branching fractions of B 0 (s) → µ + µ − µ + µ − .The combined systematic uncertainties are calculated by adding the individual components in quadrature.