Measurement of the $B^0_s\to\mu^+\mu^-$ branching fraction and effective lifetime and search for $B^0\to\mu^+\mu^-$ decays

TL;DR

This paper reports the first observation of the rare decay $B^0_s o\mu^+\mu^-$ with a measured branching fraction and lifetime, and sets an upper limit on $B^0 o\mu^+\mu^-$, all consistent with the Standard Model.

Contribution

The paper presents the first observation of $B^0_s o\mu^+\mu^-$ decay in a single experiment, including its branching fraction and effective lifetime, and provides a new upper limit for $B^0 o\mu^+\mu^-$.

Findings

Observation of $B^0_s o\mu^+\mu^-$ with 7.8 sigma significance

Measured branching fraction ${ m B}(B^0_s o\mu^+\mu^-) = (3.0 imes 10^{-9})$

First measurement of the $B^0_s o\mu^+\mu^-$ effective lifetime

Abstract

A search for the rare decays $B^0_s\to\mu^+\mu^-$ and $B^0\to\mu^+\mu^-$ is performed at the LHCb experiment using data collected in $pp$ collisions corresponding to a total integrated luminosity of 4.4 fb$^{-1}$. An excess of $B^0_s\to\mu^+\mu^-$ decays is observed with a significance of 7.8 standard deviations, representing the first observation of this decay in a single experiment. The branching fraction is measured to be ${\cal B}(B^0_s\to\mu^+\mu^-)=\left(3.0\pm 0.6^{+0.3}_{-0.2}\right)\times 10^{-9}$, where the first uncertainty is statistical and the second systematic. The first measurement of the $B^0_s\to\mu^+\mu^-$ effective lifetime, $\tau(B^0_s\to\mu^+\mu^-)=2.04\pm 0.44\pm 0.05$ ps, is reported. No significant excess of $B^0\to\mu^+\mu^-$ decays is found and a 95 % confidence level upper limit, ${\cal B}(B^0\to\mu^+\mu^-)<3.4\times 10^{-10}$, is determined. All results are in agreement with the Standard Model expectations.

Figures

Equations

$${\tau }_{{\mu }^{+}{\mu }^{-}}=\frac{{\tau }_{B_s^0}}{1-y_s^2}\left[\frac{1+2A_{\Delta \Gamma }^{{\mu }^{+}{\mu }^{-}}{y}_s+y_s^2}{1+A_{\Delta \Gamma }^{{\mu }^{+}{\mu }^{-}}{y}_s}\right],$$

$$\mathcal{B}(B_{(s)}^0\to {\mu }^{+}{\mu }^{-})=\frac{{\mathcal{B}}_{\mathrm{norm}}{ϵ}_{\mathrm{norm}}{f}_{\mathrm{norm}}}{N_{\mathrm{norm}}{ϵ}_{\mathrm{sig}}{f}_{d(s)}}\times N_{B_{(s)}^0\to {\mu }^{+}{\mu }^{-}}\equiv {\alpha }_{B_{(s)}^0\to {\mu }^{+}{\mu }^{-}}^{\mathrm{norm}}\times N_{B_{(s)}^0\to {\mu }^{+}{\mu }^{-}},$$

Full text

EUROPEAN ORGANIZATION FOR NUCLEAR RESEARCH (CERN)

​​​

CERN-EP-2017-041

LHCb-PAPER-2017-001

March 16, 2017

Measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\rightarrow\mu^{+}\mu^{-}$ decays

The LHCb collaboration†††Authors are listed at the end of this paper.

A search for the rare decays $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$ is performed at the LHCb experiment using data collected in $pp$ collisions corresponding to a total integrated luminosity of 4.4$\mbox{\,fb}^{-1}$. An excess of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decays is observed with a significance of $7.8$ standard deviations, representing the first observation of this decay in a single experiment. The branching fraction is measured to be ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})=\left(3.0\pm 0.6^{\,+0.3}_{\,-0.2}\right)\times 10^{-9}$, where the first uncertainty is statistical and the second systematic. The first measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime, $\tau({{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}})=2.04\pm 0.44\pm 0.05{\rm\,ps}$, is reported. No significant excess of $B^{0}\rightarrow\mu^{+}\mu^{-}$ decays is found and a 95% confidence level upper limit, ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})<3.4\times 10^{-10}$, is determined. All results are in agreement with the Standard Model expectations.

Published in Phys. Rev. Lett. 118 (2017), 191801

© CERN on behalf of the LHCb collaboration, licence CC-BY-4.0.

Within the Standard Model (SM) of particle physics, the ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays are very rare, because they only occur through loop diagrams and are helicity-suppressed. Since they are characterised by a purely leptonic final state, and thanks to the progress in lattice QCD calculations [1, 2, 3], their time-integrated branching fractions, ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ = $(3.65\pm 0.23)\times 10^{-9}$ and ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ = $(1.06\pm 0.09)\times 10^{-10}$ [4], are predicted in the SM with small uncertainty. These features make the $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ decays sensitive probes for physics beyond the SM, for example an extended Higgs sector [5, 6, 7]. The measurement of these processes has attracted considerable theoretical and experimental interest, culminating in the recent observation of the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay and evidence of the ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay reported by the LHCb and CMS collaborations [8]. This has been obtained by combining their datasets collected in $pp$ collisions in 2011 and 2012 [9, 10]. The measured branching fractions, ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ = $(2.8^{+0.7}_{-0.6})\times 10^{-9}$ and ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ = $(3.9^{+1.6}_{-1.4})\times 10^{-10}$, are consistent with SM predictions. The ATLAS collaboration has also recently reported a search for these decays [11].

In the ${{B}^{0}_{s}}-{{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}$ system, the light and heavy mass eigenstates are characterised by a sizable difference between their decay widths, $\Delta\Gamma=0.082\pm 0.007{\rm\,ps}^{-1}$ [12]. In the SM, only the heavy state decays to $\mu^{+}\mu^{-}$, but this condition does not necessarily hold in New Physics scenarios [13]. The contributions from the two states can be disentangled by measuring the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ effective lifetime, which, in the search for physics beyond the SM, is a complementary probe to the branching fraction measurement. The effective lifetime is defined as $\tau_{\mu^{+}\mu^{-}}\equiv\int^{\infty}_{0}t\,\Gamma\!\left(B_{s}(t)\rightarrow\mu^{+}\mu^{-}\right)dt/\int^{\infty}_{0}\Gamma\!\left(B_{s}(t)\rightarrow\mu^{+}\mu^{-}\right)dt$, where $t$ is the decay time of the ${B}^{0}_{s}$ or ${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}$ meson and $\Gamma\!\left(B_{s}(t)\rightarrow\mu^{+}\mu^{-}\right)\equiv\Gamma\!\left({{B}^{0}_{s}}(t)\rightarrow\mu^{+}\mu^{-}\right)+\Gamma\!\left({{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}(t)\rightarrow\mu^{+}\mu^{-}\right)$. The relation [14]

[TABLE]

holds, where $\tau_{{{B}^{0}_{s}}}=1.510\pm 0.005{\rm\,ps}$ is the ${B}^{0}_{s}$ mean lifetime and $y_{s}\equiv\tau_{{{B}^{0}_{s}}}\Delta\Gamma/2=0.062\pm 0.006$ [12, 15]. The parameter $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}$ is defined as $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}=-2\Re(\lambda)/(1+|\lambda|^{2})$, with $\lambda=(q/p)(A({{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}\rightarrow\mu^{+}\mu^{-})/A({{B}^{0}_{s}}\rightarrow\mu^{+}\mu^{-}))$. The complex coefficients $p$ and $q$ define the mass eigenstates of the ${{B}^{0}_{s}}-{{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}$ system in terms of the flavour eigenstates (see, e.g., Ref. [12]), and $A({{B}^{0}_{s}}\rightarrow\mu^{+}\mu^{-})$ ($A({{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}\rightarrow\mu^{+}\mu^{-}))$ is the ${{B}^{0}_{s}}$ (${{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}}$) decay amplitude. In the SM the quantity $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}$ is equal to unity but can assume any value in the range $[-1,1]$ in New Physics scenarios.

This Letter reports measurements of the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ time-integrated branching fractions, which supersede the previous LHCb results [9], and the first measurement of the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ effective lifetime. Results are based on data collected with the LHCb detector, corresponding to an integrated luminosity of 1$\mbox{\,fb}^{-1}$ of $pp$ collisions at a centre-of-mass energy $\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$, 2$\mbox{\,fb}^{-1}$ at $\sqrt{s}=8\mathrm{\,Te\kern-1.00006ptV}$ and 1.4$\mbox{\,fb}^{-1}$ recorded at $\sqrt{s}=13\mathrm{\,Te\kern-1.00006ptV}$. The first two datasets are referred to as Run 1 and the latter as Run 2.

At various stages of the analysis multivariate classifiers are employed to select the signal. In particular, after trigger and loose selection requirements, $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ candidates are classified according to their dimuon mass and the output variable, BDT, of a multivariate classifier based on a boosted decision tree [16], which is employed to separate signal and combinatorial background. The signal yield is determined from a fit to the dimuon mass distribution of candidates and is converted into a branching fraction using as normalisation modes the decays $B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{+}\rightarrow J/{\psi}K^{+}$, with $J/{\psi}\rightarrow\mu^{+}\mu^{-}$ (inclusion of charge-conjugated processes is implied throughout this Letter).

The analysis strategy is similar to that employed in Ref. [9] and has been optimised to enhance the sensitivity to both ${B}^{0}_{s}$ and ${B}^{0}$ decays to $\mu^{+}\mu^{-}$. This is achieved through a better rejection of misidentified $b$-hadron decays such as $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ (where $h^{(\prime)}=\pi,K$) and the development of an improved boosted decision tree for the BDT classifier. The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime is measured from the background-subtracted decay-time distribution of signal candidates in the lowest-background BDT region as defined later. To avoid potential biases, candidates in the dimuon mass signal region ($[5200,5445]{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$) were not examined until the analysis procedure was finalised.

The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range $2<\eta<5$, described in detail in Refs. [17, 18]. It 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 with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. Particle identification is provided by two ring-imaging Cherenkov detectors, an electromagnetic and a hadronic calorimeter, and a muon system composed of alternating layers of iron and multiwire proportional chambers. The simulated events used in this analysis are produced using the software described in Refs. [19, *Sjostrand:2007gs, *LHCb-PROC-2010-056, *Lange:2001uf, *Agostinelli:2002hh, *Allison:2006ve, *LHCb-PROC-2011-006, 26].

Candidate events for signal and normalisation are selected by a hardware trigger followed by a software trigger [27]. The $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ candidates are predominantly selected by single-muon and dimuon triggers. The $B^{+}\rightarrow J/{\psi}K^{+}$ candidates are selected in a very similar way, the only difference being a different dimuon mass requirement in the software trigger. Candidate $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ decays are used as control and normalisation channels.

The $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ candidates are reconstructed by combining two oppositely charged particles with transverse momentum with respect to the beam, $p_{\rm T}$, satisfying $0.25<p_{\rm T}<40{\mathrm{\,Ge\kern-1.00006ptV\!/}c}$, momentum $p<500{\mathrm{\,Ge\kern-1.00006ptV\!/}c}$, and high-quality muon identification [28]. Compared to the previous analysis, the muon identification requirements are tightened such that the misidentified $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ background is reduced by approximately 50%, while the signal efficiency decreases by about 10%. The muon candidates are required to form a secondary vertex with a vertex-fit $\chi^{2}$ per degree of freedom smaller than 9 and separated from any primary $pp$ interaction vertex (PV) by a flight distance significance greater than 15. Only muon candidate tracks with $\chi^{2}_{\rm IP}>25$ for any PV are selected, where $\chi^{2}_{\rm IP}$ is defined as the difference between the vertex-fit $\chi^{2}$ of the PV formed with and without the particle in question. In the selection, $B^{0}_{(s)}$ candidates must have a decay time less than $9\,\tau_{{{B}^{0}_{s}}}$, $\chi^{2}_{\rm IP}<25$ with respect to the PV for which the $\chi^{2}_{\rm IP}$ is minimal (henceforth called the $B^{0}_{(s)}$ PV), $p_{\rm T}>0.5{\mathrm{\,Ge\kern-1.00006ptV\!/}c}$ and a dimuon mass in the range $[4900,6000]{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$. A $B^{0}_{(s)}$ candidate is rejected if either of the two candidate muons combined with any other oppositely charged muon candidate in the event has a mass within 30${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ of the $J/\psi$ mass [15]. The normalisation channels are selected with almost identical requirements to those applied to the signal sample. The $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ selection is the same as that of $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$, except that the muon identification criteria are replaced with hadron identification requirements. The $B^{+}\rightarrow J/{\psi}K^{+}$ decay is reconstructed by combining a muon pair, consistent with a $J/\psi$ from a detached vertex, and a kaon candidate with $\chi^{2}_{\rm IP}>25$ for all PVs in the event. These selection criteria are completed by a loose requirement on the response of a multivariate classifier, described in Ref. [29] and unchanged since then, applied to candidates in both signal and normalisation channels. The classifier takes as input quantities related to the direction of the $B^{0}_{(s)}$ candidate, its impact parameter with respect to the $B^{0}_{(s)}$ PV, the separation between the final-state tracks, and their impact parameters with respect to any PV. After the trigger and selection requirements 78 241 signal candidates are found, which form the dataset for the subsequent branching fraction measurement.

The separation between signal and combinatorial background is achieved by means of the BDT variable, where the boosted decision tree is optimised using simulated samples of ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ events for signal and of $b\bar{b}\rightarrow\mu^{+}\mu^{-}X$ events for background. The classifier combines information from the following input variables: $\sqrt{\Delta\phi^{2}+\Delta\eta^{2}}$, where $\Delta\phi$ and $\Delta\eta$ are the azimuthal angle and pseudorapidity differences between the two muon candidates; the minimum $\chi^{2}_{\rm IP}$ of the two muons with respect to the $B^{0}_{(s)}$ PV; the angle between the $B^{0}_{(s)}$ candidate momentum and the vector joining the $B^{0}_{(s)}$ decay vertex and $B^{0}_{(s)}$ PV; the $B^{0}_{(s)}$ candidate vertex-fit $\chi^{2}$ and impact parameter significance with respect to the $B^{0}_{(s)}$ PV. In addition, two isolation variables are included, to quantify the compatibility of the other tracks in the event with originating from the same hadron decay as the signal muon candidates. Most of the combinatorial background is composed of muons originating from semileptonic $b$-hadron decays, in which other charged particles may be produced and reconstructed. The isolation variables are constructed to recognise these particles and differ in the type of tracks being considered: the first considers tracks that have been reconstructed both before and after the magnet, while the second considers tracks reconstructed only in the vertex detector. The isolation variables are determined based on the proximity of the two muon candidates to the tracks of the event and are optimised using simulated ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and $b\bar{b}\rightarrow\mu^{+}\mu^{-}X$ events. The proximity of each muon candidate to a track is measured using a multivariate classifier that takes as input quantities such as the angular and spatial separation between the muon candidate and the track, the signed distance between the muon-track vertex and the $B^{0}_{(s)}$ candidate or primary vertex, and the kinematic and impact parameter information of the track.

The BDT variable is constructed to be distributed uniformly in the range [0,1] for signal, and to peak strongly at zero for background. Its correlation with the dimuon mass is below 5%. Compared to the multivariate classifier used in the previous measurement [9], the combinatorial background with ${\rm BDT}>0.25$ is reduced by approximately 50%, mainly due to the improved performance of the isolation variables.

The expected $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ BDT distributions are determined from those of $B^{0}\rightarrow K^{+}\pi^{-}$ decays in data after correcting them for distortions due to trigger and muon identification. An additional correction is made for the ${B}^{0}_{s}$ signal, assuming the SM prediction, to account for the difference between the ${B}^{0}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ lifetimes, which affects the BDT distribution. The mass distribution of the signal decays is described by a Crystal Ball function [30]. The peak values for the ${B}^{0}_{s}$ and ${B}^{0}$ mesons are obtained from the mass distributions of $B^{0}_{s}\rightarrow K^{+}K^{-}$ and $B^{0}\rightarrow K^{+}\pi^{-}$ samples, respectively. The mass resolutions as function of $\mu^{+}\mu^{-}$ mass are determined with a power-law interpolation between the measured resolutions of charmonium and bottomonium resonances decaying into two muons. The Crystal Ball radiative tail is obtained from simulated $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ events [26], which are smeared such that they reproduce the 23${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ mass resolution measured in data.

The signal branching fractions are measured with

[TABLE]

where $N_{B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}}$ is the number of observed signal decays, $N_{\rm norm}$ is the number of normalisation-channel decays ($B^{+}\rightarrow J/{\psi}K^{+}$ and $B^{0}\rightarrow K^{+}\pi^{-}$), ${\cal B}_{\rm norm}$ is the corresponding branching fraction [15], and ${\rm\epsilon_{sig}}$ (${\rm\epsilon_{\rm norm}}$) is the total efficiency for the signal (normalisation) channel. The fraction $f_{d(s)}$ indicates the probability for a $b$ quark to fragment into a $B^{0}_{(s)}$ meson. Assuming $f_{d}=f_{u}$, the fragmentation probability $f_{\rm norm}$ for the $B^{0}$ and $B^{+}$ normalisation channel is set to $f_{d}$. The value of $f_{s}/f_{d}$ in $pp$ collision data at $\sqrt{s}=7\mathrm{\,Te\kern-1.00006ptV}$ has been measured by LHCb to be $0.259\pm 0.015$ [31]. The stability of $f_{s}/f_{d}$ at $\sqrt{s}=8\mathrm{\,Te\kern-1.00006ptV}$ and $13\mathrm{\,Te\kern-1.00006ptV}$ is evaluated by comparing the observed variation of the ratio of the efficiency-corrected yields of $B^{0}_{s}\rightarrow J/\psi\phi$ and $B^{+}\rightarrow J/{\psi}K^{+}$ decays. The effect of increased collision energy is found to be negligible for data at $\sqrt{s}=8\mathrm{\,Te\kern-1.00006ptV}$ while a scaling factor of $1.068\pm 0.046$ is applied for data at $\sqrt{s}=13\mathrm{\,Te\kern-1.00006ptV}$.

The efficiency ${\rm\epsilon_{sig(norm)}}$ includes the detector acceptance, trigger, reconstruction and selection efficiencies of the final-state particles. The acceptance, reconstruction and selection efficiencies are computed with samples of simulated events whose decay-time distributions are generated according to the SM prediction. The tracking and particle identification efficiencies are determined using control channels in data [32, 33]. The trigger efficiencies are evaluated with data-driven techniques [34].

The numbers of $B^{+}\rightarrow J/{\psi}K^{+}$ and $B^{0}\rightarrow K^{+}\pi^{-}$ decays are $(1964.2\pm 1.5)\times 10^{3}$ and $(31.3\pm 0.4)\times 10^{3}$, respectively. The normalisation factors derived from the two channels are consistent. Taking correlations into account, their weighted averages are $\alpha^{\rm norm}_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}=(5.7\pm 0.4)\times 10^{-11}$ and $\alpha^{\rm norm}_{B^{0}\rightarrow\mu^{+}\mu^{-}}=(1.60\pm 0.04)\times 10^{-11}$. In the SM scenario, the analysed data sample is expected to contain an average of $62\pm 6$ ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and $6.7\pm 0.6$ ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays in the full BDT range.

The combinatorial background is distributed almost uniformly over the mass range. In addition, the signal region and the low-mass sideband ($[4900,5200]{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$) are populated by backgrounds from exclusive $b$-hadron decays, which can be classified in two categories. The first includes $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$, $B^{0}\rightarrow\pi^{-}\mu^{+}\nu_{\mu}$, $B^{0}_{s}\rightarrow K^{-}\mu^{+}\nu_{\mu}$, and ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow p\mu^{-}\bar{\nu}_{\mu}$ decays, where one or two hadrons are misidentified as a muon. The $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$, $B^{0}\rightarrow\pi^{-}\mu^{+}\nu_{\mu}$ and ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow p\mu^{-}\bar{\nu}_{\mu}$ branching fractions are taken from Refs. [15, 35], while a theoretical estimate for $B^{0}_{s}\rightarrow K^{-}\mu^{+}\nu_{\mu}$ is obtained from Refs. [36, 37]. The mass and BDT distributions of these decays are determined from simulated samples after calibrating the $K\rightarrow\mu$, $\pi\rightarrow\mu$ and $p\rightarrow\mu$ momentum-dependent misidentification probabilities using control channels in data. An independent estimate of the $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$, $B^{0}\rightarrow\pi^{-}\mu^{+}\nu_{\mu}$ and $B^{0}_{s}\rightarrow K^{-}\mu^{+}\nu_{\mu}$ background yields is obtained by fitting the mass spectrum of $\pi^{+}\mu^{-}$ or $K^{+}\mu^{-}$ combinations selected in data, and rescaling the yields according to the $\pi\rightarrow\mu$ or $K\rightarrow\mu$ misidentification probability. The difference with respect to the results from the first method is assigned as a systematic uncertainty. The second category includes the decays $B^{+}_{c}\rightarrow J/\psi\mu^{+}\nu_{\mu}$, with $J/\psi\rightarrow\mu^{+}\mu^{-}$, and $B^{0(+)}\rightarrow\pi^{0(+)}\mu^{+}\mu^{-}$, which have at least two muons in the final state. The rate of $B^{+}_{c}\rightarrow J/\psi\mu^{+}\nu_{\mu}$ decays is evaluated from Refs. [38, 39], while those of $B^{0(+)}\rightarrow\pi^{0(+)}\mu^{+}\mu^{-}$ decays are obtained from Refs. [40, 41]. The expected yields of all exclusive backgrounds are estimated using the decay $B^{+}\rightarrow J/{\psi}K^{+}$ as the normalisation channel, with the exception of the $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ decays, which are normalised to the mode $B^{0}\rightarrow K^{+}\pi^{-}$. The contributions from $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\nu_{\mu}\bar{\nu}_{\mu}$ decays [42, 4, 43] have a negligible impact on the signal yield determination. The expected background yields with ${\rm BDT}>0.5$ in the signal region are $2.9\pm 0.3$ $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$, $1.2\pm 0.2$ $B^{+}_{c}\rightarrow J/\psi\mu^{+}\nu_{\mu}$, $0.7\pm 0.2$ ${{\mathchar 28931\relax}^{0}_{b}}\rightarrow p\mu^{-}\bar{\nu}_{\mu}$ and $0.80\pm 0.06$ $B^{0}_{(s)}\rightarrow h^{-}\mu^{+}\nu_{\mu}$ decays. The $B^{0(+)}\rightarrow\pi^{0(+)}\mu^{+}\mu^{-}$ background is negligible. Except for the misidentified $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ decays, which populate the ${B}^{0}$ signal region, the other modes are mostly concentrated in the low-mass sideband.

The Run 1 and Run 2 datasets are each divided into five subsets based on bins in the BDT variable with boundaries $0.0$, $0.25$, $0.4$, $0.5$, $0.6$ and $1.0$. The ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fractions are determined with a simultaneous unbinned maximum likelihood fit to the dimuon mass distribution in each BDT bin of the two datasets. The ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ fractional yields in each BDT bin and the parameters of the Crystal Ball functions that describe the shapes of the mass distributions are Gaussian-constrained according to their expected values and uncertainties. The combinatorial background in each BDT bin is parameterised with an exponential function, with a common slope parameter for all bins of a given dataset, while the yield is allowed to vary independently. The exclusive backgrounds are included as separate components in the fit. Their overall yields as well as the fractions in each BDT bin are Gaussian-constrained according to their expected values. Their mass shapes are determined from a simulation for each BDT bin.

The values of the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ branching fractions obtained from the fit are ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})=\left(3.0\pm 0.6^{\,+0.3}_{\,-0.2}\right)\times 10^{-9}$ and ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})=\left(1.5^{\,+1.2\,+0.2}_{\,-1.0\,-0.1}\right)\times 10^{-10}$. The statistical uncertainty is derived by repeating the fit after fixing all the fit parameters, except the $B^{0}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fractions, the background yields and the slope of the combinatorial background, to their expected values. The systematic uncertainties of ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ and ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ are dominated by the uncertainty on $f_{s}/f_{d}$ and the knowledge of the exclusive backgrounds, respectively. The correlation between the two branching fractions is negligible. The mass distribution of the $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ candidates with ${\rm BDT}>0.5$ is shown in Fig. 1, together with the fit result [44].

An excess of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ candidates with respect to the expectation from background is observed with a significance of $7.8$ standard deviations ($\sigma$), while the significance of the $B^{0}\rightarrow\mu^{+}\mu^{-}$ signal is $1.6\,\sigma$. The significances are determined, using Wilks’ theorem [45], from the difference in likelihood between fits with and without the signal component.

Since no significant ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ signal is observed, an upper limit on the branching fraction is set using the $\textrm{CL}_{\textrm{s}}$ method [46]. The ratio between the likelihoods in two hypotheses, signal plus background and background only, is used as the test statistic. The likelihoods are computed with nuisance parameters fixed to their nominal values. Pseudo-experiments are used for the evaluation of the test statistic in which the nuisance parameters are floated according to their uncertainties. The resulting upper limit on ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ is $3.4\times 10^{-10}$ at 95% confidence level.

The selection efficiency and BDT distribution of ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays depend on the lifetime, which in turn depends on the model assumption entering Eq. 1. This introduces a further model-dependence in the measured time-integrated branching fraction. In the fit, the SM value $\tau(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})=\tau_{B_{s}}/(1-y_{s})$ is assumed, corresponding to $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}=1$. The model dependence is evaluated by repeating the fit under the $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}=0$ and $-1$ hypotheses, finding an increase of the branching fraction with respect to the SM assumption of 4.6% and 10.9%, respectively. The dependence is approximately linear in the physically allowed $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}$ range.

For the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ lifetime determination, the data are background-subtracted with the sPlot technique [47], using a fit to the dimuon mass distribution to disentangle signal and background components statistically. Subsequently, a fit to the signal decay-time distribution is made with an exponential function multiplied by the acceptance function of the detector. The ${B}^{0}_{s}$ candidates are selected using criteria similar to those applied in the branching fraction analysis, the main differences being a reduced dimuon mass window, $[5320,6000]{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$, and looser particle identification requirements on the muon candidates. The former change allows the fit model for the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ signal to be simplified by removing most of the ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and exclusive background decays that populate the lower dimuon mass region, while the latter increases the signal selection efficiency. Furthermore, instead of performing a fit in bins of BDT, a requirement of BDT $>0.55$ is imposed. All these changes minimise the statistical uncertainty on the measured effective lifetime. This selection results in a final sample of 42 candidates.

The mass fit includes the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ and combinatorial background components. The parameterisations of the mass shapes are the same as used in the branching fraction analysis. The correlation between the mass and the reconstructed decay time of the selected candidates is less than 3%.

The variation of the trigger and selection efficiency with decay time is corrected for in the fit by introducing an acceptance function, determined from simulated signal events that are weighted to match the properties of the events seen in data. The use of simulated events to determine the decay-time acceptance function is validated by measuring the effective lifetime of $B^{0}\rightarrow K^{+}\pi^{-}$ decays selected in data. The measured effective lifetime is $1.52\pm 0.03{\rm\,ps}$, where the uncertainty is statistical only, consistent with the world average [15]. The statistical uncertainty on the measured $B^{0}\rightarrow K^{+}\pi^{-}$ lifetime is taken as the systematic uncertainty associated with the use of simulated events to determine the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ acceptance function.

The accuracy of the fit for the ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ effective lifetime is estimated using a large number of simulated experiments with signal and background contributions equal, on average, to those observed in the data. The contamination from ${{B}^{0}}\!\rightarrow{\mu^{+}}{\mu^{-}}$, $B^{0}_{(s)}\rightarrow h^{+}{h}^{\prime-}$ and semileptonic decays above 5320${\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ is small and not included in the fit. The effect on the effective lifetime from the unequal production rate of ${B}^{0}_{s}$ and ${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}_{s}$ mesons [48] is negligible. A bias may also arise if $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}\neq\pm 1$, with the consequence that the underlying decay-time distribution is the sum of two exponential distributions with the lifetimes of the light and heavy mass eigenstates. In this case, as the selection efficiency varies with the decay time, the returned value of the lifetime from the fit is not exactly equal to the definition of the effective lifetime even if the decay-time acceptance function is correctly accounted for. This effect has been evaluated for the scenario where there are equal contributions from both eigenstates to the decay. The result can also be biased if the background has a much longer mean lifetime than ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decays; this is mitigated by an upper decay-time cut of 13.5${\rm\,ps}$. Any remaining bias is evaluated using the background decay-time distribution of the much larger ${{B}^{0}}\!\rightarrow{{K}^{+}}{{\pi}^{-}}$ data sample. All of these effects are found to be small compared to the statistical uncertainty and combine to give 0.05${\rm\,ps}$, with the main contributions arising from the fit accuracy and the decay-time acceptance (0.03${\rm\,ps}$ each).

The mass distribution of the selected ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ candidates is shown in Fig. 2 (top). Figure 2 (bottom) shows the background-subtracted ${{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}}$ decay-time distribution with the fit function superimposed [44]. The fit results in $\tau({{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}})=2.04\pm 0.44\pm 0.05{\rm\,ps}$, where the first uncertainty is statistical and the second systematic. This measurement is consistent with the $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}=1~{}(-1)$ hypothesis at the $1.0~{}\sigma$ ($1.4~{}\sigma$) level. Although the current experimental uncertainty allows only a weak constraint to be set on the value of the $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}$ parameter in the physically allowed region, this result establishes the potential of the effective lifetime measurement in constraining New Physics scenarios with the datasets that LHCb is expected to collect in the coming years [49].

In summary, a search for the rare decays $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$ is performed in $pp$ collision data corresponding to a total integrated luminosity of 4.4$\mbox{\,fb}^{-1}$. The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal is seen with a significance of $7.8$ standard deviations and provides the first observation of this decay from a single experiment. The time-integrated $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction is measured to be $\left(3.0\pm 0.6^{\,+0.3}_{\,-0.2}\right)\times 10^{-9}$, under the $A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}=1$ hypothesis. This is the most precise measurement of this quantity to date. In addition, the first measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime, $\tau({{B}^{0}_{s}}\!\rightarrow{\mu^{+}}{\mu^{-}})=2.04\pm 0.44\pm 0.05{\rm\,ps}$, is presented. No evidence for a $B^{0}\rightarrow\mu^{+}\mu^{-}$ signal is found, and the upper limit ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})<3.4\times 10^{-10}$ at 95% confidence level is set. The results are in agreement with the SM predictions and tighten the existing constraints on possible New Physics contributions to these decays.

Acknowledgements

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); MOST and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG (Germany); INFN (Italy); NWO (The Netherlands); MNiSW and NCN (Poland); MEN/IFA (Romania); MinES and FASO (Russia); MinECo (Spain); SNSF and SER (Switzerland); NASU (Ukraine); STFC (United Kingdom); NSF (USA). We acknowledge the computing resources that are provided by CERN, IN2P3 (France), KIT and DESY (Germany), INFN (Italy), SURF (The Netherlands), PIC (Spain), GridPP (United Kingdom), RRCKI and Yandex LLC (Russia), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), PL-GRID (Poland) and OSC (USA). We are indebted to the communities behind the multiple open source software packages on which we depend. Individual groups or members have received support from AvH Foundation (Germany), EPLANET, Marie Skłodowska-Curie Actions and ERC (European Union), Conseil Général de Haute-Savoie, Labex ENIGMASS and OCEVU, Région Auvergne (France), RFBR and Yandex LLC (Russia), GVA, XuntaGal and GENCAT (Spain), Herchel Smith Fund, The Royal Society, Royal Commission for the Exhibition of 1851 and the Leverhulme Trust (United Kingdom).