###### Abstract

This article reviews the most recent measurements of
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decay properties at the LHC, which are
the most precise to date. The measurements of the branching fraction and
effective lifetime of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay by the
ATLAS, CMS and LHCb collaborations, as well as the search for
$B^{0}\rightarrow\mu^{+}\mu^{-}$ decays are summarised with a focus on the
experimental challenges. Furthermore prospects are given for these
measurements and new observables that become accessible with the foreseen
amounts of data by the end of the LHC.

###### keywords:

$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$, FCNC, Rare Decay, Flavour physics,
ATLAS, CMS, LHCb, LHC

1 1 0 https://doi.org/ Analysis of $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$
decays at the LHC Title Kai-Feng Chen 1,†*, Titus Mombächer 2,†* and Umberto
de Sanctis 3,†* Kai-Feng Chen, Titus Mombächer and Umberto de Sanctis Chen,
K.-F.; Mombächer, T.; de Sanctis, U. Correspondence: KC: Kai-
<EMAIL_ADDRESS>TM<EMAIL_ADDRESS>US:
<EMAIL_ADDRESS>The authors contributed equally to this work.

##  1 Introduction

This review summarises the most recent measurements related to the
$B^{0}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$
decays performed with the ATLAS, CMS and LHCb experiments.

The $B^{0}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$
decays belong to the category of Flavour Changing Neutral Current (FCNC)
processes and therefore highly suppressed in the Standard Model (SM). This
makes them important tools in the search for New Physics (NP) since they can
provide indirect constraints on NP processes that interfere with the SM
processes and alter the rates and decay properties sizeably. They are even
sensitive to particles that are out of the kinematic range accessible by
particle colliders, including the Large Hadron Collider (LHC). The
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays are among the most sensitive
FCNC processes due to their small theoretical uncertainty and clean
experimental signature Bobeth et al. (2014a, b); Hermann et al. (2013); Beneke
et al. (2018, 2019). In the SM, the decays
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ are forbidden at leading-order and can
proceed only via loop diagrams. In addition, they are also suppressed by the
helicity conservation and the presence of off-diagonal CKM matrix elements,
leading to very small expected decay time integrated branching fractions.
Additional interest in $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decays comes from
the simple description in effective field theory Altmannshofer et al. (2012);
Beaujean et al. (2012). The decays can only proceed via axial-vector (Wilson
coefficient $\mathcal{C}_{10}$), scalar ($\mathcal{C}_{S}$) or pseudo-scalar
($\mathcal{C}_{P}$) $b\rightarrow s\ell\bar{\ell}$ currents, where the scalar
and pseudo-scalar currents are forbidden in the SM. Thus, measurements of
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ properties are crucial inputs for global
fits of the parameters that govern $b\rightarrow s\ell\bar{\ell}$ transitions.

The most up-to-date SM predictions for the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$
branching fractions are calculated in Ref. Beneke et al. (2019) and yield

$\begin{array}[]{lll}{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})&=&(3.66\pm
0.14)\times 10^{-9}\,\quad\textrm{and}\\\ {\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})&=&(1.03\pm 0.05)\times
10^{-10}\,\quad.\end{array}$ (1)

The values of these branching fractions are $C\\!P$ averaged and time-
integrated and include final state radiation effects, so that they can be
readily compared with the experimental measurements, which do not distinguish
between $B^{0}_{(s)}$ $C\\!P$ eigenstates and take final state radiation
effects into account. Next-to-leading order electroweak corrections and next-
to-next-to-leading order QCD corrections are also included in the
calculations. Recently, several progresses in lattice quantum chromodynamics
(QCD) Aoki et al. (2019); Bazavov et al. (2018); Bussone et al. (2016);
Dowdall et al. (2013); Hughes et al. (2018), in the calculation of electroweak
effects at next-to-leading order Bobeth et al. (2014b), and QCD effects at
next-to-next-to-leading order Hermann et al. (2013) helped significantly in
reducing the theoretical uncertainties on both branching fractions. Enhanced
electromagnetic contributions from virtual photon exchange were also proven to
produce larger corrections to theoretical uncertainties than previously
assumed Beneke et al. (2018, 2019). These predictions also take into account
the finite width difference measured in the $B^{0}_{s}$ system, that apply in
the experimental measurements where the data samples are un-tagged (see Ref.
De Bruyn et al. (2012); Bruyn et al. (2012)). Alternative predictions are also
available. They are obtained using the relation between
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays and $\Delta m_{d(s)}$, the mass
difference of the $B^{0}_{(s)}$ mass eigenstates Buras (2003); King et al.
(2019). In addition, it has been recently pointed out Buras (2023) that the
current way to calculate ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$
could be affected by the presence on NP effects. Therefore, a calculation
based on $\Delta m_{d(s)}$ and $|\epsilon_{K}|$ considering only the SM
contribution has been proposed. In both cases, the resulting values for ${\cal
B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ are slightly different than the ones
shown in Eq. 1, but still compatible within the theoretical uncertainties of
the calculation, that are, on their own, still smaller than the experimental
precision. Since all Collaborations used the values reported in Eq. 1 to
assess the level of compatibility of the measurements with the SM predictions,
in the remainder of the article, the values quoted in Eq. 1 will be used as
reference values for the ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ SM
predictions.

While the mentioned reference does not quote a value for the ratio of the two
branching fractions, this can be easily calculated as:

${\mathcal{R}}=\frac{{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})}{{\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})}=\frac{\tau_{B^{0}}}{1/\Gamma_{H}^{s}}\left(\frac{f_{B^{0}}}{f_{B^{0}_{s}}}\right)^{2}\left|\frac{V_{td}}{V_{ts}}\right|^{2}\tfrac{M_{B^{0}}\sqrt{1-\tfrac{4m_{\mu}^{2}}{M_{B^{0}}^{2}}}}{M_{B^{0}_{s}}\sqrt{1-\tfrac{4m_{\mu}^{2}}{M_{B^{0}_{s}}^{2}}}}=0.0281\pm
0.0016\quad$ (2)

where $\tau_{B^{0}}$ and $1/\Gamma_{H}^{s}$ are the lifetimes of the $B^{0}$
and of the heavy mass eigenstate of the $B^{0}_{s}$; $M_{B^{0}_{s}}$ and
$M_{B^{0}}$ are the masses and $f_{B^{0}_{s}}$ and $f_{B^{0}}$ the meson decay
constants of the $B^{0}_{s}$ and $B^{0}$ mesons respectively; $V_{td}$ and
$V_{ts}$ the elements of the CKM matrix and $m_{\mu}$ the mass of the muon.
Using the same input values as Ref.Beneke et al. (2019), the numerical value
in Eq. 2 is obtained. It is worth noting that the ratio has a theoretical
uncertainty which is smaller than the single branching fractions due to the
cancellation of most of the factors. In particular this ratio has the same
value in all theories obeying the Minimal Flavour Violation (MFV) paradigm
(including the SM) and as such it is a test of the latter. It is therefore of
additional interest to evaluate ${\mathcal{R}}$ also in upcoming measurements.

A second observable of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay
considered in the latest experimental results is its effective lifetime
$\tau_{\mu\mu}$. This observable is complementary to the branching fraction
because it is sensitive to potential New Physics (NP) effects which are
flavour-dependent De Bruyn et al. (2012). In fact, in the SM, only the heavy
$C\\!P$ odd heavy-mass eigenstate component of the
$B^{0}_{s}$-$\overline{B}^{0}_{s}$ system contributes to the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay amplitude: an assumption which does
not hold in every NP scenario. Therefore, the measurement of this quantity
could reveal the presence of NP effects which do not affect the branching
fraction measurement. $\tau_{\mu\mu}$ is simply defined as the mean lifetime
of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decays

$\displaystyle\tau_{\mu\mu}$ $\displaystyle\equiv$
$\displaystyle\frac{\int^{\infty}_{0}t\left<\Gamma\left(B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\right)\right>dt}{\int^{\infty}_{0}\left<\Gamma\left(B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\right)\right>dt}$
(3) $\displaystyle=$
$\displaystyle\frac{\tau_{B^{0}_{s}}}{1-y_{s}^{2}}\left[\frac{1+2\mathcal{A}_{\Delta\Gamma}y_{s}+y_{s}^{2}}{1+\mathcal{A}_{\Delta\Gamma}y_{s}}\right],$

where $t$ is the proper decay time of the $B^{0}_{s}$ meson and $y_{s}$ and
the $C\\!P$ parameter $\mathcal{A}_{\Delta\Gamma}$ are defined as

$y_{s}\equiv\frac{\Delta\Gamma_{s}}{2\Gamma_{s}}\qquad\mathcal{A}_{\Delta\Gamma}\equiv\frac{R_{H}^{\mu^{+}\mu^{-}}-R_{L}^{\mu^{+}\mu^{-}}}{R_{H}^{\mu^{+}\mu^{-}}+R_{L}^{\mu^{+}\mu^{-}}}$
(4)

and $R_{H}^{\mu^{+}\mu^{-}}$ and $R_{L}^{\mu^{+}\mu^{-}}$ denote the
contributions of the heavy and light mass eigenstates of the $B^{0}_{s}$
system to the un-tagged $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay rate. Since
the $\mu^{+}\mu^{-}$ final state is $C\\!P$-odd, in the SM
$\mathcal{A}_{\Delta\Gamma}=+1$ and the effective lifetime is equal to the
lifetime of the heavy-mass $B^{0}_{s}$ eigenstate $\tau^{SM}_{\mu\mu}$. The
$C\\!P$ asymmetry $\mathcal{A}_{\Delta\Gamma}$ can receive contributions from
NP effects, particularly from scalar and pseudoscalar operators, even in cases
where the branching fraction ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ is
not modified. The most recent $\tau^{SM}_{\mu\mu}$ value is $1.624\pm
0.009{\rm\,ps}$ Amhis et al. (2023), which can be slightly different from that
used by the various Collaborations depending on publication time of their most
recent measurement.

All experimental results described in this review assume the SM hypothesis
$\mathcal{A}_{\Delta\Gamma}=1$ in the calculation of efficiencies and
acceptance for the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay and thus for its
branching fractions. However, the ATLAS, CMS and LHCb Collaborations estimated
the impact of such assumption on the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$
branching fraction, which spans from 4% to 10% depending on
$\mathcal{A}_{\Delta\Gamma}$ varying in the interval [-1,1].

The article is structured as follows: in Sections 2, 3 and 4, the measurements
by the ATLAS, CMS and LHCb Collaborations are reported respectively, while in
Section 5 the results obtained by the latest official LHC combination are
presented. Section 6 provides a summary of the status of the measurements and
prospects of the three Collaborations for the HL-LHC phase.

##  2 The ATLAS ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ and effective
lifetime measurements

###  2.1 The ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ measurement

The measurement of the branching fractions of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$
decays performed by the ATLAS Collaboration is described and documented in
Ref. ATLAS Collaboration (2019). The analysis uses 26.1 $\mbox{\,fb}^{-1}$ of
Run-2 data collected at $\sqrt{s}$ = 13 $\mathrm{\,Te\kern-1.00006ptV}$, and
combines the result with the previously published Run-1 analysis ATLAS
Collaboration (2016) on 4.7 $\mbox{\,fb}^{-1}$ of data at $\sqrt{s}$ = 7
$\mathrm{\,Te\kern-1.00006ptV}$ and 20.3 $\mbox{\,fb}^{-1}$ at $\sqrt{s}$ = 8
$\mathrm{\,Te\kern-1.00006ptV}$.

In order to remove the dependence from the knowledge of the $b$-quark
production cross section and minimise the systematic uncertainties, the
branching fractions are measured relative to a reference channel. For its
abundance and well-measured branching ratio, the $B^{+}\rightarrow
J/\psi(\mu^{+}\mu^{-})K^{+}$ decay channel has been chosen for this purpose.
As a consequence, the procedure to extract the ${\cal
B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ takes into account the difference in
the fragmentation fractions $f_{u,d,s}$ of $b$-quarks to form, respectively, a
$B^{+}$, $B^{0}_{s}$ or $B^{0}$ meson. Also the different acceptances and
efficiencies between the signal and the reference channels are taken into
account. Hence, the branching fractions ${\cal
B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ are expressed as:

${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})={N_{B^{0}_{(s)}}\over
N_{B^{+}}}{f_{u}\over f_{d(s)}}{\epsilon_{\rm tot}^{B^{+}}\over\epsilon_{\rm
tot}}{\cal B}(B^{+}\rightarrow J/\psi(\mu^{+}\mu^{-})K^{+})~{},$ (5)

where $N_{B^{0}_{(s)}}$ ($N_{B^{+}}$) is the measured yield of
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ ($B^{+}\rightarrow J/\psi K^{+}$)
events, and $\epsilon_{\rm tot}$ ($\epsilon_{\rm tot}^{B^{+}}$) is the total
signal ($B^{+}\rightarrow J/\psi K^{+}$) efficiency. Events from
$B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\phi$
decay, with $J/\psi\rightarrow\mu\mu$ and $\phi\rightarrow KK$, are also used
as control sample for the signal kinematic variables exploited in the
analysis.

The signal selection starts with a hardware dimuon trigger requiring one muon
with transverse momentum $p_{\rm T}$ $>$ 4 $\mathrm{\,Ge\kern-1.00006ptV}$ and
the other with $p_{\rm T}$ $>$ 6 $\mathrm{\,Ge\kern-1.00006ptV}$. In the
offline analysis, both muons are required to have the same $p_{\rm T}$
thresholds as in the trigger selection, to have pseudo-rapidity
$|\eta|~{}<~{}$2.5, and to pass stringent track-quality requirements (Tight
muons). Signal candidates are formed with two muons with opposite electric
charges. Kaon candidates for the reference channel are reconstructed in the
tracking system and are required to have $p_{\rm T}$ $>$ 1
$\mathrm{\,Ge\kern-1.00006ptV}$ and $|\eta|~{}<~{}$2.5.

$B$-meson kinematic observables are reconstructed imposing quality
requirements on the dimuon vertex for the signal, or on the vertex formed by
the dimuon system and one track for the reference channel. The reconstructed
$B$ candidates are also required to fall within a fiducial volume defined as
$\mbox{$p_{\rm T}$}(B)$ $>$ 8 $\mathrm{\,Ge\kern-1.00006ptV}$ and
$|\eta(B)|~{}<~{}$2.5.

The analysis uses mainly the $B$-candidate invariant mass to characterise the
selected events. $B$-candidates with a mass in the 4766-5966 MeV interval, are
considered. A blind analysis is performed where the dimuon invariant mass
signal region between 5166 and 5526 $\mathrm{\,Me\kern-1.00006ptV}$ is not
used until the analysis criteria and strategies are finalised.

The main backgrounds for this analysis can be split into three categories:
continuum background, partially reconstructed $B$ decays (PRD) and peaking
background. The continuum background consists mainly of muons produced in
uncorrelated hadron decays. It is the dominant background for the analysis and
it is several orders of magnitude larger than the signal. Therefore, a Boosted
Decision Tree Hoecker et al. (2007) (c-BDT) is emplyed to efficiently separate
the signal from this background type. The c-BDT is based on 15 kinematic
variables with high discriminating power which describe the kinematics of the
$B$-meson candidate, the secondary vertex displacement, the kinematic
properties of the muons and the rest of the event (such as the isolation of
the $B$ candidate and that of the two muon tracks with respect to the rest of
the event). The c-BDT is trained and validated on the data mass sidebands.

The PRD background is made of decays where the two muons in the final state
come from one of the following topologies: $(a)$ ’cascade’ transitions with
the muons coming from the same ancestor (e.g. $b\rightarrow c\mu\nu\rightarrow
s\mu\mu\nu\nu$), and labelled same-side muons (SS); $(b)$ from the same decay
chain (e.g. $B\rightarrow J/\psi X$ or $B^{0}\rightarrow\mu\mu K^{*}$) and
labelled same vertex muons (SV); $(c)$ from $B_{c}\rightarrow J/\psi\mu\nu$
decays; $(d)$ from semileptonic $B$ decays where a hadron $h$ ($\pi$, $K$ or
proton) is misidentified as a muon (e.g. $B\rightarrow\mu h\nu$). All these
types of backgrounds populate the low-mass sideband with contributions also
into the dimuon mass signal region.

The peaking background consists of charmless two body decays
$B^{0}_{(s)}\rightarrow h^{+}h^{(^{\prime})-}$ ($h^{(\prime)}$ being a pion or
a kaon) that are reconstructed as signal events due to the hadrons being
misidentified as muons. This background component falls in the signal region
and presents the same features of the $B^{0}\rightarrow\mu^{+}\mu^{-}$ signal.
Its contribution has been studied with the help of a dedicated MC sample and
validated in data in a region enriched by hadrons misidentified as muons. The
resulting peaking background contribution is estimated to be $2.9\pm 2.0$
events in the signal region.

To extract the ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ using Eq. 5,
the yield of the reference channel and the efficiency ratio between the two
channels needs to be computed. The
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ yield
$N_{B^{+}}$ is obtained by an unbinned extended maximum-likelihood fit to the
$\mu\mu K^{+}$ invariant mass distribution, where the shape parameters are
fitted simultaneously in data and simulation.

The efficiency ratio between signal and reference channels is computed from
appropriate simulation samples within the fiducial volume of the analysis.
These samples are reweighted in such a way that they reproduce the
distributions of the number of primary vertices (and therefore pile-up),
$\mbox{$p_{\rm T}$}(B)$, $|\eta(B)|$ and trigger efficiencies (as a function
of $\mbox{$p_{\rm T}$}(\mu)$ and $|\eta(\mu)|$) as measured in data.
Furthermore, a correction to the $B^{0}_{s}$ lifetime in the simulated signal
sample is applied to match the distribution of the heavy $B^{0}_{s}$ mass
eigenstate, because the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay proceeds in
the SM exclusively through the heavy $B^{0}_{s}$ mass eigenstate, as described
in Section 1.

The yields of signal events $N_{B^{0}_{(s)}}$ are extracted simultaneously
from an unbinned extended maximum-likelihood fit to the dimuon invariant mass
distribution $m_{\mu\mu}$. In order to enhance the sensitivity of the
analysis, four bins in the c-BDT output (in increasing order of signal-over-
background ratio) are defined in order to have constant signal efficiency
equal to 18% in each bin. The fit is performed simultaneously in the four
c-BDT bins. The first c-BDT bin, which has the lowest signal-over-background
ratio, is dominated by the main backgrounds. It is introduced in the fit to
improve the backgrounds modelling and reduce the systematic uncertainties
related to them. The $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ signals are
parameterised by a double Gaussian function to take into account different
resolutions in the dimuon invariant mass depending on the different regions of
the ATLAS detector. The shape and the relative signal efficiencies are assumed
to be the same in all c-BDT bins. The continuum background is described by a
first order polynomial, while the background coming from SS and SV events is
parameterised with an exponential function. These backgrounds are fluctuated
independently in each c-BDT bin. Finally, the description of the peaking
background is based on the same model used to describe the signal, with a
constraint on the total yield of $2.9\pm 2.0$ equally distributed in the c-BDT
bins.

The ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ values are extracted
through a simultaneous unbinned extended maximum-likelihood fit using the
components written in Eq. 5 and the $N_{B^{0}_{(s)}}$ event yields extracted
from the invariant mass fits just described. The ${\cal B}(B^{+}\rightarrow
J/\psi(\mu^{+}\mu^{-})K^{+})$ value is taken as the world average from the PDG
Tanabashi et al. (2018), while the hadronisation probability ratio
$\frac{f_{u}}{f_{d(s)}}=0.256\pm 0.013$ from the HFLAV average Amhis et al.
(2017).

The measurements are dominated by statistical uncertainties, with the most
prominent sources of systematic uncertainty coming from: the fit uncertainties
(where the largest contributors are the mass scale and the
$b\rightarrow\mu^{+}\mu^{-}X$ background parameterisation), the
$\frac{f_{u}}{f_{d(s)}}$ ratio (only for the BR($B^{0}_{s}$) measurement) and
the reference channel yield. All systematic uncertainties are described in the
likelihood as Gaussian constraints.

A Neyman construction Neyman  is employed to extract the 68.3%, 95.5% and
99.7% confidence intervals in the ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$–${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ plane. The likelihood function from the
described Run-2 result is then combined with the likelihood function from the
Run-1 result ATLAS Collaboration (2016). The only common parameters in the
combination are: the fitted ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$
and the external inputs (${\cal B}(B^{+}\rightarrow
J/\psi(\mu^{+}\mu^{-})K^{+})$ and $\frac{f_{u}}{f_{d(s)}}$). All remaining
nuisance parameters are treated as uncorrelated between the two results.

The ATLAS results, obtained by combining 25 $\mbox{\,fb}^{-1}$ from Run1 and
26.1 $\mbox{\,fb}^{-1}$ from Run2 LHC campaigns, are ATLAS Collaboration
(2019):

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ $\displaystyle=$
$\displaystyle\left(2.8\,^{+0.8}_{-0.7}\right)\times 10^{-9}\quad,$ (6)
$\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ $\displaystyle=$
$\displaystyle\left(-1.9\pm 1.6\right)\times 10^{-10}\quad,$

with a significance for the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal of 4.6
standard deviations ($\sigma$). The 95 % confidence level (CL) upper limit for
the $B^{0}\rightarrow\mu^{+}\mu^{-}$ signal is ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})<2.1\times 10^{-10}$, as obtained with the
Neyman procedure described in Ref. Neyman . Fig. 1 shows the dimuon invariant
mass distribution in the highest-score BDT bin (left) and the likelihood
contours in the ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ \- ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ plane (right).

Figure 1: (Left): Dimuon invariant mass distributions in data, for the bin
with the highest scores of BDT output. Superimposed is the result of the
maximum-likelihood fit. The total fit is shown as a continuous line, with the
dashed lines corresponding to the observed signal component, the
$b\rightarrow\mu^{+}\mu^{-}X$ background, and the continuum background. The
signal components are grouped in one single curve, including both the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and the (negative)
$B^{0}\rightarrow\mu^{+}\mu^{-}$ component. The curve representing the peaking
$B^{0}_{(s)}\rightarrow hh^{\prime}$ background lies very close to the
horizontal axis. (Right): Likelihood contours for the combination of the Run 1
and 2015-2016 Run 2 results (shaded areas). The contours are obtained from the
combined likelihoods of the two analyses, for values of
$-2\Delta\mathrm{ln}(L)$ equal to 2.3, 6.2 and 11.8. The empty contours
represent the result from 2015-2016 Run 2 data alone. The SM prediction with
uncertainties is also indicated. Figures from Ref. ATLAS Collaboration (2019).

###  2.2 The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime
measurement

Using the same dataset and the same configurations for the event selection and
the simulated samples, ATLAS has subsequently performed a measurement of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime $\tau_{\mu\mu}$ ATLAS
Collaboration (2023). As explained in Section 1, the measurement of this
quantity is complementary to the branching fraction measurement in the
searches for NP phenomena. The only difference between the two analyses lies
in the different selection applied to the c-BDT output. A requirement on the
c-BDT output to be larger than 0.365 is applied to the dataset, while all
other requirements are the same as the BR analysis. The value of this
requirement was selected after an optimisation procedure based on the
maximisation of the $S/\sqrt{S+B}$ figure-of-merit.

The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime is measured using
a binned $\chi^{2}$ fit to the proper decay time distribution of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal component after the subtraction of
the background. The proper decay time $\tilde{t}_{\mu^{+}\mu^{-}}$ is defined
as
$\tilde{t}_{\mu^{+}\mu^{-}}=\frac{L_{xy}m_{B_{s}^{0}}^{\mathrm{PDG}}}{p_{\mathrm{T}}^{B_{s}^{0}}}$,
where $L_{xy}$ is the decay length projected along the reconstructed
$B_{s}^{0}$ momentum in the transverse plane, $m_{B_{s}^{0}}^{\mathrm{PDG}}$
the world averaged mass of $B_{s}^{0}$ mesons from Ref. Workman and Others
(2022) and $p_{\mathrm{T}}^{B_{s}^{0}}$ the magnitude of the candidate’s
reconstructed transverse momentum. To extract $\tau_{\mu\mu}$, three main
steps have been completed:

  * •

A fit to the dimuon invariant mass, in the same range as for the BR analysis

  * •

The extraction of the $\tilde{t}_{\mu^{+}\mu^{-}}$ distribution of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ component using the sPlot technique Pivk
and Le Diberder (2005)

  * •

A binned $\chi^{2}$ fit to $\tilde{t}_{\mu^{+}\mu^{-}}$ distribution comparing
Monte-Carlo simulated effective lifetime templates corresponding to different
values of $\tau_{\mu\mu}$.

In the first step, the dimuon invariant mass distribution, after all selection
cuts described above, is fit using a five parameters model made of three
Probability Density Functions (PDF): a double Gaussian to describe the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ component, a linear function to describe
the combinatorial (or continuum) background component and an exponential
function to describe the PRD component. Additional resonant and non-resonant
backgrounds (such as $B\rightarrow hh^{\prime}$, $B_{c}^{\pm}$ and
semileptonic $B$ decays), as well as the $B^{0}\rightarrow\mu^{+}\mu^{-}$
component, are neglected in this fit and considered as sources of systematic
uncertainties whose impact on $\tau_{\mu\mu}$ is evaluated through MC pseudo-
experiments (as described later in the text). The fit yields 58$\pm$13 events
in the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ mass window.

The second step exploits the sPlot statistical technique to extract the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal proper-decay time component from
the invariant mass fit. The signal proper-decay time distribution is
background-subtracted by means of per-event weights computed using the result
of the invariant mass fit described above.

The third and final step consists of a binned-$\chi^{2}$ fit to the proper-
decay time distribution extracted in the previous step. This distribution is
considered in the interval 0-12 ps and divided in twelve equal width bins.
Pure signal proper-decay time simulated templates in the same interval and
binning scheme corresponding to different values of $\tau_{\mu\mu}$ are
generated, and a $\chi^{2}$-binned fit is performed with respect to
background-subtracted data. The $\chi^{2}$ calculation takes both the
statistical uncertainty on the weight-corrected MC and the Poissonian
uncertainty in each data bin as expected from the predicted MC content for
that bin into account. The template minimising the $\chi^{2}$ corresponds to
an observed lifetime $\tau_{\mu\mu}^{\mathrm{Obs}}$ of 1.07 ps. MC pseudo-
experiments studies, generated for a lifetime of 1.624 ps (i.e. the SM
predicted value) showed that the lifetime extraction procedure had a bias of
82 fs due to the low-statistics regime of the fit. This bias is found to be
constant in the $B^{0}_{s}$ lifetime range considered in the analysis.
Therefore the quoted value for $\tau_{\mu\mu}^{\mathrm{Obs}}$ has been
corrected for this effect. The statistical uncertainty on
$\tau_{\mu\mu}^{\mathrm{Obs}}$ is instead extracted using a MC pseudo-
experiments based Neyman construction, yielding to a value of
$\tau_{\mu\mu}^{\mathrm{Obs}}=0.99^{+0.42}_{-0.07}(\textrm{stat.})$ ps. Fig. 2
shows the signal proper decay time distribution extracted from data
superimposed with the MC template minimising the $\chi^{2}$ distribution
(left) and the MC pseudo-experiments based Neyman construction used to
estimate the statistical uncertainty of the measurement (right).

The dominant systematic uncertainties for this measurement are related to the
data-MC discrepancies (134 fs evaluated in data by repeating, under the same
statistical regime as the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal case,
the same fit procedure in the $B^{\pm}\rightarrow J/\psi K^{\pm}$ channel), to
the background mass and lifetime models (86 fs), to the fit dependence from
the lifetime used in MC pseudo-experiments generation and the $B^{0}_{s}$
eigenstates admixture (15 fs evaluated by generating MC pseudo-experiments in
the $\tau^{SM}_{L}$ \- $\tau^{SM}_{H}$ lifetime interval) and the neglected
resonant and non-resonant backgrounds (12 fs). The total systematic
uncertainty is then obtained by summing in quadrature and symmetrising the
impact on $\tau_{\mu\mu}$ of all single sources. This yields to an observed
value of $\tau_{\mu\mu}^{\mathrm{Obs}}$ of
$0.99^{+0.42}_{-0.07}(\textrm{stat.})\pm 0.17(\textrm{syst.})$. The value is
compatible with the SM prediction of 1.624 ps ($\mathcal{A}_{\Delta\Gamma}=1$)
as well as with the other experimental results described in this article.

Figure 2: (Left): Signal proper decay time distribution extracted with the
sPlot background subtraction procedure applied to the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ invariant mass fit. The superimposed
signal MC template is the result of the lifetime fit procedure discussed in
the text. The uncertainties on the data points are calculated as Poisson
fluctuations on the MC yield prediction (continuous red histogram) in the
corresponding bin. (Right): 68% and 95% CL bands obtained with a Neyman
construction based on MC pseudo-experiments for the signal and background
components. The yellow lines interpolate the band boundaries in order to
smooth the effects of limited number of MC pseudo-experiments used. The
dashed-dotted blue line corresponds to the average expected
$\tau_{\mu\mu}^{\mathrm{Obs}}$ value at a given $\tau_{\mu\mu}$value used to
generate MC pseudo-experiments. The horizontal dashed black line corresponds
to the experimentally observed value of $\tau_{\mu\mu}^{\mathrm{Obs}}$=0.99
ps, yielding a 68% CL band of [0.92,1.41] ps (thick vertical dashed purple
lines) and a 95% CL band of [0.77,1.73] ps (thin vertical dashed purple
lines). The same construction at the
$\tau_{\mu\mu}^{\mathrm{Obs}}$corresponding to $\tau_{\mu\mu}$= 1.624 ps (the
SM prediction) yields [1.44,2.26] ps as 68% CL band. Figures from Ref. ATLAS
Collaboration (2023).

##  3 Measurement of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay properties and
search for $B^{0}\rightarrow\mu^{+}\mu^{-}$ decay at CMS

The latest analysis by the CMS collaboration is based on the LHC Run-2 data
collected in 2016–2018 at a center-of-mass energy of 13 TeV, corresponding to
an integrated luminosity of 140 fb-1 CMS Collaboration (2022). The studies
based on LHC Run-1 samples collected in 2011–2012 can be found in the earlier
publications CMS Collaboration (2020). There is no attempt to combine the
latest publication with the results from 2011–2012 data as the expected gain
in sensitivity is modest. In this section the latest CMS measurement of ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$, the search for
$B^{0}\rightarrow\mu^{+}\mu^{-}$ decay, and the effective lifetime measurement
using $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ events are discussed.

###  3.1 Measurement of ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ and
search for $B^{0}\rightarrow\mu^{+}\mu^{-}$ decay

The characteristic signal $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ comprises two
muons originating from a single displaced vertex, isolated from other
activities, with momentum aligned with the flight direction, and an invariant
mass peaking at $M(B^{0}_{s})$ or $M(B^{0})$. The primary contributors to the
background comprise combinatorial events, involving instances where the two
muons originate from different heavy quarks, partially reconstructed
semileptonic decays wherein both muons emanate from the same $B$ meson (with
one of the muons from a misidentified charged hadron), and the background
arising from peaking charmless two-body hadronic $B$ meson decays.

The data events were collected with a set of dimuon triggers for this study:
the L1 trigger required two oppositely charged muons within the range of
$|\eta|<1.5$, while at HLT the dimuon should form a secondary vertex and are
required to be within specific mass ranges. The dimuon candidates are used to
reconstruct $B$ mesons for the signal and normalization
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ and
$B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\phi$
channels. The selections are reserved to be as similar as possible for partial
cancellation of systematic effects. Muons at offline analysis are required to
have a high-quality track fit at tracker, a transverse momentum at least 4 GeV
and $|\eta|<1.4$. To suppress misidentified muons from charged pion and kaon
decays, a multivariate-analysis(MVA)-based algorithm has been introduced.
Extra kaons are required in the reconstruction for the normalization channels.
A trajectory representing the $B$ candidate is built from the decay vertex and
$B$ candidate’s momentum, and is extrapolated to the closest point for each
reconstructed primary vertex; the primary vertex with the smallest distance to
the extrapolated point is selected for the analysis.

How to reduce the combinatorial and partially reconstructed backgrounds are
the main challenges to the study. To enhance the analysis sensitivity, a
dedicated MVA discriminator, combining various discriminating observables into
a single score distribution ($d_{\rm MVA}$) using a boosted decision tree
algorithm, is introduced. The inputs for $d_{\rm MVA}$ includes pointing
angles, defined as the angles between the B momentum and the direction
connecting the primary and secondary vertices, observables related to the
secondary vertex such as quality of the vertex finding, and observables that
are designed to identify nearby decay products in semi-leptonic decays of $b$
and $c$ hadrons. The $d_{\rm MVA}$ training is employed by an advanced
gradient boosting algorithm, supported by the XGBoost library Chen and
Guestrin (2016). The training utilises a mix of
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal events and background events
selected from the data sidebands. Following a fine-tuning of input observables
to align the kinematics of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ decays
(considering variations in the uncertainties of the dimuon vertex position),
the control decay $B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}K^{+}$ channel is employed to evaluate the performance of the $d_{\rm
MVA}$ in data.

Charmless two-body decays $B^{0}_{(s)}\rightarrow h^{+}h^{(^{\prime})-}$, like
$B^{0}\rightarrow K^{+}\pi^{-}$ and $B^{0}_{s}\rightarrow K^{+}K^{-}$, can
mimic the signal when both charged hadrons are misidentified as muons. The
misidentification probabilities in data are calculated by utilizing
$K^{0}_{S}\rightarrow\pi^{+}\pi^{-}$, $\phi(1020)\rightarrow K^{+}K^{-}$, and
$\Lambda\rightarrow p\pi^{-}$ decays, restricting the decay distance of
$K^{0}_{S}$ and $\Lambda$ candidates to align with the lifetime of the $B$
meson. Misidentification of pions and kaons primarily originates from their
decays into muons. An agreement between the observed data and simulations is
observed for both pions and kaons. The proton misidentifying rate is much
smaller hence the contributions from the associated processes are totally
negligible. After stringent multivariate-based muon identification
requirement, the charmless two-body backgrounds reduce to a negligible level.

Because of the limited precision in measuring the b-quark production cross
section at the LHC, directly determining the branching fraction (${\cal
B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$) could introduce significant
uncertainty. As a common practice, the signal branching fraction is assessed
by normalizing it to the
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ decay
channel. In addition the
$B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\phi$
decays, with $J/\psi\rightarrow\mu\mu$ and $\phi\rightarrow KK$, are
considered as a cross-check, and might become more precise if the ${\cal
B}(B^{0}_{s}\rightarrow J/\psi\phi(1020))$ is further improved by future
B-factory studies. Another advantage of measuring branching fractions in a
relative manner is the potential cancellation of systematic uncertainties
common in the selection of the signal and normalization channels. The exact
formulae for the $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ branching fractions
are similar to those used in Eq.5:

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle={\cal B}(B^{+}\rightarrow J/\psi
K^{+}){N_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}\over N_{B^{+}\rightarrow J/\psi
K^{+}}}{\epsilon_{B^{+}\rightarrow J/\psi
K^{+}}\over\epsilon_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}}{f_{u}\over f_{s}},$
(7) $\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle={\cal
B}(B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\phi){N_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}\over
N_{B_{s}^{0}\rightarrow J/\psi\phi}}{\epsilon_{B^{+}\rightarrow
J/\psi\phi}\over\epsilon_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}},$ (8)
$\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ $\displaystyle={\cal
B}(B^{+}\rightarrow J/\psi K^{+}){N_{B^{0}\rightarrow\mu^{+}\mu^{-}}\over
N_{B^{+}\rightarrow J/\psi K^{+}}}{\epsilon_{B^{+}\rightarrow J/\psi
K^{+}}\over\epsilon_{B^{0}\rightarrow\mu^{+}\mu^{-}}}{f_{u}\over f_{d}},$ (9)

where the yields and the selection efficiencies for each processes are denoted
by $N_{X}$ and $\epsilon_{X}$ ($X=B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$,
$B^{0}\rightarrow\mu^{+}\mu^{-}$, $B^{+}\rightarrow J/\psi K^{+}$, or
$B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\phi$).
The production fractions for $B^{+}$, $B^{0}$, and $B^{0}_{s}$ mesons are
represented by $f_{u}$, $f_{d}$, and $f_{s}$. The ratio $\frac{f_{u}}{f_{d}}$
is set to unity due to isospin symmetry, while the ratio
$\frac{f_{u}}{f_{s}}$, together with ${\cal B}(B^{+}\rightarrow J/\psi K^{+})$
and ${\cal B}(B^{0}_{s}\rightarrow J/\psi\phi)$, are external inputs.

The results are obtained with a simultaneous unbinned maximum likelihood fits
in multiple categories. For the measurement of branching fractions, a two-
dimensional fit using the dimuon invariant mass and its uncertainty as
observables is introduced. The events are categorized according to data-taking
period, signal purity based on $d_{\rm MVA}$, and $|\eta|$ of the most-forward
muon. The likelihood function include five components:
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$
signals, semileptonic background, peaking two-body decays, and the
combinatorial events. The signal components are represented using Crystal Ball
functions for the dimuon mass. The width of these Crystal Ball functions is
parameterized based on the per-event mass resolution. To model the mass
resolution, a kernel estimation approach is employed, utilizing Gaussian
kernels. The semileptonic background is modeled by a Gaussian with free
parameters in the fit to the data, while the peaking background is modeled by
a sum of Gaussian and Crystal Ball functions with the shape parameters
determined from simulated events. The yields of the semileptonic and peaking
background components are first derived and then included in the fit with
uncertainties from the hadron to muon misidentifying rate as constrained
nuisance parameters. The combinatorial background is modeled by a linear
function with yields and slope free to vary in the fits.

For the branching fraction measurements the experimental uncertainties include
signal efficiency corrections due to mismodeling of $d_{\rm MVA}$, the charged
kaon efficiency in the normalization channels, trigger efficiencies, and
fitting bias, while the rest of uncertainties are smaller than 1%. The
mismodeling of the $d_{\rm MVA}$ distribution has been investigated through
two distinct studies with $B^{+}\rightarrow J/\psi K^{+}$ events. In the first
study, a direct comparison is conducted between background-subtracted data,
with the sPlot technique Pivk and Le Diberder (2005) on the $B^{+}\rightarrow
J/\psi K^{+}$ invariant mass distribution, and the simulated distributions.
The second study involves reweighting of the simulated samples to align with
the data distributions, employing the XGBoost tool. The disparity between the
two studies is quantified as a systematic uncertainty. The systematic
uncertainty arising from the selection of background models is derived through
pseudo-experiments, incorporating variations in the fit. The uncertainties in
the input branching fractions of the normalization channels and the
$\frac{f_{u}}{f_{d(s)}}$ ratio, are implemented as external uncertainties.

The resulting branching fractions for $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}\rightarrow\mu^{+}\mu^{-}$ are:

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[3.83^{+0.38}_{-0.36}\,{\rm(stat)}^{+0.19}_{-0.16}\,{\rm(syst)}^{+0.14}_{-0.13}\,{(f_{s}/f_{u})}\right]\times
10^{-9},$ (10) $\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[0.37^{+0.75}_{-0.67}\,{\rm(stat)}^{+0.08}_{-0.09}\,{\rm(syst)}\right]\times
10^{-10}.$ (11)

The results incorporate external inputs, specifically ${\cal
B}(B^{+}\rightarrow J/\psi K^{+})=(1.020\pm 0.019)\times 10^{-3}$, ${\cal
B}(J/\psi\rightarrow\mu^{+}\mu^{-})=(5.961\pm 0.033)\times 10^{-2}$, and
$f_{s}/f_{u}=0.231\pm 0.008$. The input $f_{s}/f_{u}$ value is derived from
the $p_{T}$-dependent measurement by LHCb Aaij et al. (2021) and the $p_{T}$
distribution observed in this measurement. Figure 3 shows the dimuon invariant
mass distributions from the categories with different signal purity; the
results of the fit are superimposed. The profile likelihood contours enclose
the regions with different coverage are shown in Figure 4. Alternatively the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction is measured using the
$B^{0}_{s}\rightarrow J/\psi\phi$ decays as the normalization, which leads to

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[4.02^{+0.40}_{-0.38}\,{\rm(stat)}^{+0.28}_{-0.23}\,{\rm(syst)}^{+0.18}_{-0.15}\,{(\mathcal{B})}\right]\times
10^{-9},$ (12)

where the last uncertainty arises from the uncertainty in the
$B^{0}_{s}\rightarrow J/\psi\phi$ branching fraction (${\cal
B}(B^{0}_{s}\rightarrow J/\psi\phi)=(1.04\pm 0.040)\times 10^{-3}$). The
lifetime of the $B^{0}_{s}$ meson has a significant impact on the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction too; a scaling factor
on the resulting branching fraction ($1.577-0.358\cdot\tau_{B^{0}_{s}}$, where
$\tau_{B^{0}_{s}}$ is $B^{0}_{s}$ lifetime in ps) for alternative lifetime
hypotheses other than the SM value (1.61 ps) is provided. The upper limit on
the $B^{0}\rightarrow\mu^{+}\mu^{-}$ decay is calculated to be ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})<1.9\times 10^{-10}$ at 95% confidence
level, using the CLs method Read (2002).

Figure 3: The dimuon invariant mass distributions for the candidates with
$d_{\rm MVA}>0.99$ (left) and $0.99>d_{\rm MVA}>0.90$ (right) categories. The
solid blue curves are the projections of fit model, while the individual
components of the fit are also presented. Figures from Ref. CMS Collaboration
(2022).  Figure 4: The profile likelihood (left) in ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ versus ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ 2D plane, while the contours enclose the
regions with $1\sigma$-$5\sigma$ coverage. Figure from Ref. CMS Collaboration
(2022).

###  3.2 Measurement of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ Effective
Lifetime

The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime $\tau_{\mu\mu}$is
extracted with a unbinned maximum likelihood in three-dimensions including
dimuon invariant mass, decay time, and decay time uncertainty. The decay time
$t_{\mu^{+}\mu^{-}}$, which is calculated for each event, is defined by the
product of the flight length and the invariant mass of the $B$ candidate,
divided by the magnitude of the $B$ candidate momentum. The events are also
categorized in the data-taking period, purity based on $d_{\rm MVA}$, and the
pseudorapidity of the most forward muon. The dimuon invariant mass
distribution is modeled with the same functions introduced for the branching
fraction measurements, while the decay time distribution for signal events is
modeled by an exponential function convoluted with the decay time resolution
function. The decay time resolution function is parameterized with the
measured decay time uncertainty. The acceptance as a function of the decay
time is obtained from simulated events and corrected with the
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ events
from data. The decay time distribution for combinatorial background decays is
obtained from high-mass sideband events. The decay time uncertainty models
used in the fit are obtained from simulation samples and mass sideband data as
well.

The systematic uncertainties in the lifetime measurement are mostly driven by
the correlations between the $d_{\rm MVA}$ and the decay time, as the key
input variables for the $d_{\rm MVA}$ classifier: the pointing angle of $B$
candidate and its associated uncertainty are strongly correlated with the
decay time observable. Any mismodeling in the simulation results in
significant impacts on the decay time distribution. A correction as a ratio of
the decay time distributions for different $d_{\rm MVA}$ requirements is
derived from $B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}K^{+}$ events. This method introduced a bias up to 0.1 ps for the data
recorded in 2016, and reduced in the later data sets. The possible bias arises
in fitting and modeling is also tested with
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ events,
but with a relaxed selection criterion. Other systematic uncertainties are
minor, estimated to be smaller than 0.01 ps.

The resulting effective lifetime for $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$
events is:

$\displaystyle\tau_{\mu\mu}=1.83^{+0.23}_{-0.20}\,{\rm(stat)}^{+0.04}_{-0.04}\,{\rm(syst)}~{}{\rm
ps},$ (13)

which is consistent with the SM prediction and the other experimental results
described in this article. The decay time distribution for the candidates in
the region of $5.28<m_{\mu^{+}\mu^{-}}<5.48$ GeV is shown in Fig. 5.

Figure 5: The proper decay time distribution for the candidates in the region
of $5.28<m_{\mu^{+}\mu^{-}}<5.48$ GeV, with the result of the fit
superimposed. The solid blue curve is the sum of all fit component, while the
shaded areas are the background components. Figure from Ref. CMS Collaboration
(2022).

##  4 Analysis of $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays with LHCb

The most recent analysis of $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ with the
LHCb experiment Aaij et al. (2022a, b) was performed with the full
$pp$-collision data collected in the LHC Run 1 and Run 2 campaigns. The total
integrated luminosity corresponds to $1\,\text{fb}^{-1}$ at
$\sqrt{s}=7\,\mathrm{\,Te\kern-1.00006ptV}$, $2\,\text{fb}^{-1}$ at
$\sqrt{s}=8\,\mathrm{\,Te\kern-1.00006ptV}$ and $6\,\text{fb}^{-1}$ at
$\sqrt{s}=13\,\mathrm{\,Te\kern-1.00006ptV}$.

In total the analysis comprises the search and branching fraction measurements
of the decays $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$,
$B^{0}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ with initial state radiation
($B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ was only investigated in the
region $m(\mu^{+}\mu^{-})>4.9\,{\mathrm{\,Ge\kern-1.00006ptV\\!/}c^{2}}$), as
well as the measurement of the effective lifetime of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay. A precise branching fraction
measurement is achieved by normalising the signal decay with two high-
statistics decay modes, $B^{0}\rightarrow K^{+}\pi^{-}$ and
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ with
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\\!\rightarrow\mu^{+}\mu^{-}$,
similarly to what is done by ATLAS and CMS Collaborations and shown in Eq.5
for the $B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$
channel. The decay modes $B^{0}\rightarrow K^{+}\pi^{-}$ and
$B^{0}_{s}\rightarrow K^{+}K^{-}$ are used as control modes for the effective
lifetime measurement as well.

Dominant background processes mimicking the signal on the one hand arise from
random combinations of two muons from two different $b$-hadron decays in the
same event. On the other hand they can come from $b$-hadron decays where one
or more final state particles have been wrongly identified as a muon.
Furthermore, $b$-hadron decays wher part of the decay products have not been
reconstructed can constitute a background. The selection of the signal decays
largely inherits from previous analyses of a subset of the data Aaij et al.
(2017) and targets particularly the selecting of
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays over aforementioned backgrounds,
whereas the measurement of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ is a
byproduct of the analysis. The LHCb detector, as used to collect the above
mentioned data, employed a two-staged online selection. Firstly, events are
selected by a hardware trigger that requires at least one muon with a high
transverse momentum. Secondly, a two-staged software trigger is applied, which
performs a full event reconstruction. In the software trigger, events
fulfilling minimum requirements on the muon momentum and its impact parameter,
are kept. Also events are kept where these requirements are met by non-signal
candidates to maximise signal efficiency.

In the offline selection, candidate $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$
decays are selected by combining two well reconstructed oppositely charged
particles identified as muons Archilli et al. (2013) with a transverse
momentum in the range of $0.25<\mbox{$p_{\rm
T}$}<40\,{\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. The common vertex is required
to have a good vertex fit quality and be clearly separated from the associated
$pp$-collision vertex. The resulting $B^{0}_{(s)}$ candidate is required to
have a transverse momentum greater than
$0.5\,{\mathrm{\,Ge\kern-1.00006ptV\\!/}c}$. Candidates in the full
instrumented pseudorapidity region $2<\eta<5$ are retained for analysis. A
preliminary selection based on a Boosted Decision Tree (BDT) is applied to
remove a large fraction of combinatorial background while maintaining a high
signal efficiency. The BDT is trained with variables related to the decay
topology of two particles originating from a vertex displaced with respect to
the primary vertex. A highly efficient veto on the combination of a signal
muon with another particle in the event identified as muon that result in a
dimuon mass close to the ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ mass
allows to effectively remove
$B_{c}^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\mu^{+}\nu$ decays. A selection on a combination of particle
identification algorithms is performed and tuned to maximise the
$B^{0}\rightarrow\mu^{+}\mu^{-}$ significance ($h=K,\pi$), suppressing
$B^{0}_{(s)}\rightarrow h^{+}h^{-}$ and $\mathchar
28931\relax^{0}_{b}\rightarrow p\mu^{-}\nu$ decays. The final selection is
performed on a second BDT, called in the following s-BDT. This s-BDT includes,
apart from variables related to the decay topology, notably isolation
classifiers - specifically developed for this analysis - that inspect the
closeness of the signal muon tracks to other tracks in the event that are
either reconstructed in all tracking detector stations or only in the detector
closest to the collision region. The $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$
yields are measured by fitting the dimuon invariant mass distribution in bins
of this final selection s-BDT, discarding only the lowest bin (that
corresponds to about $25\,\%$ of the signal) in order to maximise the signal
sensitivity. The samples of $B^{0}\rightarrow K^{+}\pi^{-}$ and
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ are
selected in a similar way except for trigger and particle identification
criteria for the $B^{0}\rightarrow K^{+}\pi^{-}$ mode and removing the
${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$ veto. For $B^{0}\rightarrow
K^{+}\pi^{-}$, the muon identification criteria are replaced by hadron
identification and a trigger selection independent of the candidate is
required to achieve an unbiased selection.

###  4.1 Measurement of the branching fractions of
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$

In order to achieve unbiased branching fraction estimates, efficiencies are
calculated either on corrected simulation or directly on data. Importantly,
the fractions of the s-BDT bins are determined from
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ simulation, where the $B^{0}_{(s)}$
quantities and the number of tracks in the event are reweighted from data-
simulation comparisons in high-statistics
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$ and
$B^{0}_{s}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}\phi$
samples. The resulting corrected s-BDT fractions are then independently cross-
checked with $B^{0}\rightarrow K^{+}\pi^{-}$ data samples, corrected by the
different trigger and particle identification response. Measuring the
branching fractions relative to two modes, $B^{0}\rightarrow K^{+}\pi^{-}$ and
$B^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}K^{+}$, allows
for a stringent cross check of the efficiencies by calculating the ratio
between the estimated branching fractions of the two and comparing it to the
ratio of the published branching fractions Zyla et al. (2020). An excellent
agreement is found.

The invariant mass shape of signal $B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$
decays is described with two-sided Crystal Ball functions Skwarnicki (1986),
where the mean of the Gaussian core is calibrated from $B^{0}_{s}\rightarrow
K^{+}K^{-}$ and $B^{0}\rightarrow K^{+}\pi^{-}$ data samples. The mass
resolution of about $22\,{\mathrm{\,Me\kern-1.00006ptV\\!/}c^{2}}$ is
determined from the interpolation of the measured resolutions of charmonium
and bottomonium resonances. The tail parameters are estimated from simulation.
Small differences in the resolution and the tail parameters are found to
appear across the s-BDT bins and are accounted for in the final fit.

Exclusive background decays remaining in the fully selected samples have been
carefully studied with simulation, calibrated in data. A large focus in the
most recent analysis is laid on the correct estimation of the
misidentification of charged hadrons as muons. Decays of the form
$B^{0}_{(s)}\rightarrow h^{+}h^{(^{\prime})-}$ ($h=K,\pi$) with both charged
hadrons misidentified create a peaking structure very close to the
$B^{0}\rightarrow\mu^{+}\mu^{-}$ peak and therefore form the most relevant
remaining background component. Misidentification occurs in the detector
dominantly because the hadrons decay in-flight into muons. The hadron
misidentification rate is estimated with a dedicated procedure using
$D^{0}\\!\rightarrow K^{-}\pi^{+}$ from $D^{*+}\\!\rightarrow D^{0}\pi^{+}$
decays from simulation and data. This procedure takes explicitly into account
that the $D^{0}\\!\rightarrow K^{-}\pi^{+}$ invariant mass shape deforms
significantly with hadrons decaying in-flight. As additional cross check, the
misidentification rate is investigated from $B^{0}\rightarrow K^{+}\pi^{-}$
data samples by determining the $B^{0}\rightarrow K^{+}\pi^{-}$ yield in $\pi$
$\mu$, $K$ $\mu$ and $\pi$ $K$ mass distributions.

A summary of the final mass fit to obtain the signal branching fractions is
displayed in Fig. 6.

Figure 6: Left: mass distribution of the selected
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ candidates (black dots) with s-BDT
$>0.5$. The result of the fit is overlaid and the different components are
detailed: $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ (red solid line),
$B^{0}\rightarrow\mu^{+}\mu^{-}$ (green solid line),
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ (violet solid line), combinatorial
background (blue dashed line), $B^{0}_{(s)}\rightarrow h^{+}h^{-}$ (magenta
dashed line), $B^{0}\\!\rightarrow\pi^{-}\mu^{+}\nu_{\mu}$,
$B^{0}_{s}\\!\rightarrow K^{-}\mu^{+}\nu_{\mu}$,
$B_{c}^{+}\\!\rightarrow{J\mskip-3.0mu/\mskip-2.0mu\psi\mskip
2.0mu}\mu^{+}\nu_{\mu}$ and $\mathchar 28931\relax^{0}_{b}\\!\rightarrow
p\mu^{-}\nu_{\mu}$ (orange dashed line), and
$B^{0(+)}\\!\rightarrow\pi^{0(+)}\mu^{+}\mu^{-}$ (cyan dashed line). The solid
bands around the signal shapes represent the variation of the branching
fractions by their total uncertainty. Right: two-dimensional profile
likelihood of the branching fractions for the
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays. The measured central values of
the branching fractions are indicated with a blue dot. The profile likelihood
contours for $68\,\%$, $95\,\%$ and $99\,\%$ CL regions of the result are
shown as blue contours, while the brown contours indicate the previous
measurement Aaij et al. (2017) and the red cross shows the SM prediction.
Figures from Ref. Aaij et al. (2022b).

A precise measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching
fraction is obtained to be

$\displaystyle{\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})=\left(3.09^{+0.46}_{-0.43}{}^{+0.15}_{-0.11}\right)\times
10^{-9},$ (14)

where the first uncertainties are of statistical and the second of systematic
nature. The systematic uncertainties are dominated by the knowledge of the
ratio of fragmentation fractions $\frac{f_{u}}{f_{d(s)}}$ of $B^{0}_{s}$ and
$B^{0}$ mesons which enters the normalisation equation because the decay is
measured relative to $B^{0}$ and $B^{+}$ decays. The
$B^{0}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ decays are not observed and
consequently upper limits on their branching fractions are set to

$\displaystyle\mathcal{B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle<1.2\times 10^{-10}\text{, and}$ (15)
$\displaystyle\mathcal{B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$
$\displaystyle)<2.0\times 10^{-9}$ (16)

at $95\,\%$ CL, respectively. Similarly, an upper limit on the branching
fraction ratio $\mathcal{R}$ was determined at $95\,\%$ CL to

$\displaystyle\mathcal{R}<0.095$ (17)

These values include systematic uncertainties, which are dominated by the
knowledge of the background components that include misidentified hadrons. A
correlation of $11\,\%$ is observed between the measurement of the
$B^{0}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$
components.

###  4.2 Measurement of the effective lifetime of the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay

The effective lifetime of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay has
been measured on the same sample with a slightly different selection. Since
there is effectively no background from hadron-muon misidentification in the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ mass peak region, the dimuon mass window
is adapted to exclude these backgrounds and the particle identification
requirements are loosened to increase the signal yield. The conditions of
triggered events are required to be met either from the signal candidate
itself or the remainder of the event, which facilitates the modelling of the
acceptance. Furthermore the data are analysed in only two bins of the final
selection s-BDT, chosen to maximise the sensitivity to the effective lifetime.
The mass distributions in each s-BDT region are fitted independently to
extract background-subtracted decay time distributions with the sPlot
technique Pivk and Le Diberder (2005). A simultaneous fit to the two
background-subtracted decay-time distributions as shown in Fig. 7 is employed
to extract the effective lifetime. In order to extract an unbiased lifetime
measurement, the acceptance effects of the reconstruction selection
requirements have to be modelled. The decay time acceptance is modelled by
fitting parametric functions to the efficiency distribution in simulation,
where the simulation has been weighted to improve data-simulation differences.
The procedure is validated by measuring the lifetimes of $B^{0}_{s}\rightarrow
K^{+}K^{-}$ and $B^{0}\rightarrow K^{+}\pi^{-}$ in data, finding good
agreement with the world average values Zyla et al. (2020). The uncertainty of
the measurement of the $B^{0}_{s}\rightarrow K^{+}K^{-}$ lifetime is taken as
systematic uncertainty. Further systematic effects like the sample
contamination from $B^{0}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}_{(s)}\rightarrow h^{+}h^{-}$ decays, acceptance modelling,
uncertainties in the background decay time distributions and
$B^{0}_{s}$-$\overline{B}^{0}_{s}$-production asymmetries are investigated and
are found to have only sub-leading to negligible effects. The measured
effective lifetime is found to be

$\displaystyle\tau_{\mu^{+}\mu^{-}}=2.07\pm 0.29\pm 0.03\,\text{ps},$ (18)

where the first uncertainty is statistical and the second systematic. This
value is outside the lifetime interval defined by the $B^{0}_{s}$ light
($\mathcal{A}_{\Delta\Gamma}=-1$) and heavy ($\mathcal{A}_{\Delta\Gamma}=1$)
mass eigenstates, but is consistent with these values at the level of 2.2 and
1.5 standard deviations, respectively.

Figure 7: The background-subtracted decay-time distributions with the fit
model used to determine the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective
lifetime superimposed. The distributions in the low and high BDT regions are
shown in the left and right plot, respectively. Figures from Ref. Aaij et al.
(2022b)

##  5 Combination of the measurements by the LHC experiments

The latest combination of the measurements from LHC experiments is presented
in Ref. ATL (2020). It includes the results from ATLAS ATLAS Collaboration
(2019), CMS Sirunyan et al. (2019), and LHCb Aaij et al. (2017). The
combination is performed based on the two-dimensional profile likelihoods
obtained by each experiment from the fits to the dimuon events. Such a method
allows to properly take into account the correlations between ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ and ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$, and the upper limit of ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ can be evaluated using the same inputs.
Note the previous combination paper Khachatryan et al. (2015) based on CMS and
LHCb data collected during LHC Run-1 is based on an unbinned maximum
likelihood fit simultaneously to the events from both experiments.

ATLAS results, described in Section 2, are extracted from the data samples of
25 fb-1 collected at the center-of-mass energies $\sqrt{s}=$ 7 and 8 TeV, and
26.3 fb-1 at $\sqrt{s}=$ 13 TeV. The resulting branching fractions were
obtained:

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[2.8^{+0.8}_{-0.7}\right]\times 10^{-9},$ (19)
$\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[-1.9\pm 1.6\right]\times 10^{-10},$ (20)

where the systematic uncertainties are included in the evaluation. The
corresponding significance for $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ signal is
$4.6\sigma$, while the upper limit for $B^{0}\rightarrow\mu^{+}\mu^{-}$
branching fraction is $<2.1\times 10^{-10}$ at 95% CL.

CMS analysis is based on the data samples of 5 fb-1, 20 fb-1, and 36 fb-1,
collected at the center-of-mass energies of $\sqrt{s}=$ 7, 8 and 13 TeV,
respectively Sirunyan et al. (2019). The resulting branching fractions and
effective lifetime for $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ are

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[2.9^{+0.7}_{-0.6}\,{\rm(exp)}\pm
0.2\,{\rm(frag)}\right]\times 10^{-9},$ (21) $\displaystyle{\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[0.8^{+1.4}_{-1.3}\right]\times 10^{-10},$ (22)
$\displaystyle\tau_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}$
$\displaystyle=1.70^{+0.60}_{-0.43}\,{\rm(stat)}\pm 0.09\,{\rm(syst)}~{}{\rm
ps}.$ (23)

The signals for $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and
$B^{0}\rightarrow\mu^{+}\mu^{-}$ yield a significance of $5.6\sigma$ and
$1.0\sigma$, respectively. The upper limit for
$B^{0}\rightarrow\mu^{+}\mu^{-}$ branching fraction is evaluated to be
$<3.6\times 10^{-10}$ at 95% CL. The first uncertainty of ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ combined statistical and systematic
uncertainties from the analysis, while the second uncertainty is from the
uncertainty in the fragmentation ratio $f_{d}/f_{s}$.

LHCb studies are performed on the data samples of 1 fb-1, 2 fb-1, and 1.4
fb-1, collected at the center-of-mass energies $\sqrt{s}=$ 7, 8 and 13 TeV,
respectively Aaij et al. (2017). The analysis yields the following results:

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[3.0\pm
0.6\,{\rm(stat)}^{+0.3}_{-0.2}\,{\rm(syst)}\right]\times 10^{-9},$ (24)
$\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[1.5^{+1.2}_{-1.0}\,{\rm(stat)}^{+0.2}_{-0.1}\,{\rm(syst)}\right]\times
10^{-10},$ (25) $\displaystyle\tau_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}$
$\displaystyle=2.04\pm 0.44\,{\rm(stat)}\pm 0.05\,{\rm(syst)}~{}{\rm ps},$
(26)

with signal significances of $7.8\sigma$ and $1.6\sigma$ for
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}\rightarrow\mu^{+}\mu^{-}$
decays, respectively. An upper limit ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})<3.4\times 10^{-10}$ at 95% CL is obtained.

For the combination of decay branching fractions, profiled likelihoods are
computed in the two-dimensional grid of ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ and ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ plane from each experiment and the SM
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ lifetime is assumed. As the current
measurements are dominating by statistical uncertainties, the systematic
uncertainties are treated independently for the three measurements, except for
the common nuisance parameter, $f_{d}/f_{s}$ ratio. The $f_{d}/f_{s}$
uncertainty is profiled separately in each likelihood and retained only in the
LHCb experiment. To test the impact of this correlation, the ${\cal
B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ are evaluated with and without the
$f_{d}/f_{s}$ uncertainty in ATLAS and CMS liikelihoods. The impact is found
to be negligible. Additionally, the dependence of $f_{d}/f_{s}$ on the
transverse momentum is checked and is found to be consistent within the
assigned uncertainties.

The profiled likelihood for each experiment is then modeled with a two-
dimensional variable-width Gaussian, which describes asymmetric likelihoods
(and asymmetric uncertainties) and also the correlation between the two
branching fractions. This analytical function is found to be consistent with
the original likelihood for each experiment. The log-likelihoods from the
three measurements are summed across the ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ \- ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ grid points and then fitted using the
variable-width Gaussians. By maximizing the modeled likelihood function the
combined branching fractions and the associated uncertainties are derived:

$\displaystyle{\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[2.69^{+0.37}_{-0.35}\right]\times 10^{-9},$ (27)
$\displaystyle{\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$
$\displaystyle=\left[0.6\pm 0.7\right]\times 10^{-10}.$ (28)

The upper limit on ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ is evaluated as
$<1.6(1.9)\times 10^{-10}$ at 90% (95%) CL, which is calculated under the
positive ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ hypothesis by
renormalising the likelihood in the interested region. The combined ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ branching fraction is found to be
lower than any single result, which is due to the strong anti-correlation
between two branching fractions. The individual profiled likelihood (left) and
the combined likelihood in the ${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
\- ${\cal B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ plane (right) are shown in Fig.
8.

The compatibility with the SM predictions is estimated to be 2.4$\sigma$ for
${\cal B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$, 0.64$\sigma$ for ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$, and 2.1$\sigma$ if computed in the ${\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$ \- ${\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})$ plane. These values are calculated
assuming Wilks’ theorem and with theoretical uncertainties included. In
addition to the individual branching fractions, a combined estimation on the
ratio of branching fractions $\mathcal{R}$ (see Eq. 2 is also derived:

$\displaystyle\mathcal{R}$ $\displaystyle={{\cal
B}(B^{0}\rightarrow\mu^{+}\mu^{-})\over{\cal
B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})}=0.021^{+0.030}_{-0.025},$ (29)

where the corresponding upper limit is evaluated to be $\mathcal{R}<0.052$
($0.060$) at 90% (95%) CL.

Figure 8:  Left plot: the two-dimensional likelihood contours for the
$B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-}$ decays from ATLAS (red dashed line),
CMS (green dot-dashed line), and LHCb (blue long-dashed line) experiments,
together with contours for their combination (continuous line). The likelihood
contours are corresponding to the values of $-2\Delta\ln\mathcal{L}$ = 2.3,
6.2, and 11.8, respectively. Right plot: the likelihood contours for the
combination of the three results, corresponding to the values of
$-2\Delta\ln\mathcal{L}$ = 2.3, 6.2, 11.8, 19.3, and 30.2, or 1 to 5 $\sigma$
levels in a bidimensional Gaussian approximation. Figures from Ref. ATL
(2020).

The $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime is measured in the
last analysis iteration by all three experiments, as reported in Sections 2–4.
However, at the time when the combination was perforemed, only CMS and LHCb
Collaborations had a measurement of this quantity. Therefore, a combination
has been carried out based only on their results, exploiting a similar method
as for the ${\cal B}(B_{(s)}^{0}\rightarrow\mu^{+}\mu^{-})$ combination. The
LHCb analysis is carried out with a bin-likelihood fit to the background-
subtracted decay time distribution, while the CMS measurement is carried out
with a two-dimensional likelihood fit to the decay time and dimuon invariant
mass distributions. As the analyses are fully dominated by the statistical
uncertainties, the combination is performed by describing CMS and LHCb
likelihoods (as a function of effective lifetime
$\tau_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}$) with variable-width Gaussians,
and then, to determine their combined value, the two measurements are assumed
to be uncorrelated. The resulting $\tau_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}$
value and the corresponding uncertainty are:

$\displaystyle\tau_{B^{0}_{s}\rightarrow\mu^{+}\mu^{-}}$
$\displaystyle=1.91^{+0.37}_{-0.35}~{}{\rm ps}.$ (30)

Both CMS and LHCb Collaborations have recently released updated analyses, as
discussed in Sections 3-4. Another iteration of the combination is foreseen in
the near future, incorporating the results from all three experiments based on
the full Run 2 LHC campaign data.

##  6 Conclusion and prospects

In recent years the ATLAS, CMS and LHCb Collaborations made a push towards
precision measurements of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching
fraction, which resulted in measurements that reach a precision of down to
$10\,\%$ relative uncertainty. These measurements are the most precise to
date. At the same time all three Collaborations have begun measuring the
effective lifetime of the decay to understand the $C\\!P$ structure of the
decay. Contrary to initial evidence in the first combination of CMS and LHCb
measurements Khachatryan et al. (2015), the $B^{0}\rightarrow\mu^{+}\mu^{-}$
decay has not been confirmed yet. All results are in good agreement with the
SM, strongly constraining potential NP scenarios. To achieve even higher
sensitivities, a community effort is ongoing to combine the results of all
three experiments. The results of the previous combination have been presented
in this review, but have been superseded by the legacy measurements by the CMS
and LHCb Collaborations. Once the measurement with the full Run 2 data of the
ATLAS Collaboration is published as well, this combination will be repeated to
have the most precise picture possible with the harvest of Run 2 data.

After the LHC Run 2, in 2021 the experiments began to take data again with
increased instantaneous luminosity until the end of 2025. After that, the
High-Luminosity LHC phase will begin, which will have increased pile-up
conditions for all experiments and a massively increased total luminosity. The
ATLAS and CMS experiments will strongly upgrade their detectors to cope with
the increased pile-up conditions. However, they also target a significant
dimuon mass resolution improvement by $20\,\%-30\,\%$ (ATLAS) and
$40\,\%-50\,\%$ (CMS), respectively. The LHCb experiment is planning to follow
and go through a major upgrade in 2031 to begin taking data with the LHC Run
5. By the end of the LHC lifetime, ATLAS and CMS aim to have collected
$3000\,\mbox{\,fb}^{-1}$, while LHCb is estimating $300\,\mbox{\,fb}^{-1}$.
Under these conditions and assuming the central values as predicted by the
Standard Model, the ATLAS, CMS and LHCb collaborations made extrapolations to
the expected sensitivity of future measurements ATL (2022); LHC (2018). For
the ATLAS experiment, the sensitivity strongly depends on the trigger
conditions for dimuon events with the upgraded detector. In the most
conservative scenario the expected statistical-only (statistical and
systematic) uncertainties reach $19\,\%$ ($23\,\%$) relative to the central
value for the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction and
$134\,\%$ ($135\,\%$) for the $B^{0}\rightarrow\mu^{+}\mu^{-}$ branching
fraction, while in the most optimistic scenario they reach $5\,\%$ ($13$) for
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $25\,\%$ ($26\,\%$) for
$B^{0}\rightarrow\mu^{+}\mu^{-}$. The dominant systematic uncertainties in
these projections arise from external inputs like the uncertainty on the
fragmentation fraction ratio $f_{s}/f_{d}$ and the branching fraction
uncertainty of the normalisation channel.

The CMS collaboration expects to reach uncertainties of $7\,\%$ on the
branching fraction of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ and $16\,\%$ on the
branching fraction of $B^{0}\rightarrow\mu^{+}\mu^{-}$. The expected
uncertainty on the effective $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ lifetime is
$0.05\,$ps. This precision will allow stringent constraints on the parameter
$A^{\mu^{+}\mu^{-}}_{\Delta\Gamma}$ and in particular break the degeneracy
between possible scalar and pseudoscalar contributions beyond the SM to this
decay.

The LHCb collaboration expects to reach a statistical uncertainty on the
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction of $1.8\%$, however,
the analysis will be systematically limited by the external inputs of the
fragmentation fraction ratios and the normalisation branching fractions, which
are estimated to become $4\,\%$ by then. On the contrary, the ratio
$\mathcal{B}(B^{0}\rightarrow\mu^{+}\mu^{-})/\mathcal{B}(B^{0}_{s}\rightarrow\mu^{+}\mu^{-})$
is not expected to become systematically limited and is expected to reach a
relative precision of $10\,\%$. The measurement of the effective
$B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ lifetime is expected to reach a precision
of $0.033\,$ps. Both the CMS and LHCb Collaborations expect to establish the
$B^{0}\rightarrow\mu^{+}\mu^{-}$ decay signal at more than $5\sigma$ level.

The expected large yield of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decays will
also allow to access the $C\\!P$ parameter $S_{\mu^{+}\mu^{-}}$, which
describes the time-dependent $C\\!P$-violation in the decay Buras et al.
(2013). Adding this parameter will complete the base of $C\\!P$ observables
and provide complementary constraints to physics beyond the SM that are not
constrained by the other observables. A nonzero value of this parameter will
be an immediate sign for a $C\\!P$-violating phase beyond the SM. This
parameter can only be determined by measuring the decay-time distribution of
$B^{0}_{s}$ and $\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}_{s}$ decays
separately and thus requires the tagging of the $B^{0}_{s}$ flavour. Assuming
a similar performance of the flavour tagging as in Run 2, the LHCb
collaboration expects to reach a precision of $0.2$ of this parameter.
Provided a sufficient flavour tagging performance can be achieved, this
analysis could potentially be performed by the CMS and ATLAS experiments.

To achieve the projected sensitivities discussed in this section and possibly
surpass them, it will be important to maintain the basic assumptions. For the
ATLAS and CMS experiments it will be crucial to design trigger strategies that
allow to keep the muon transverse momentum thresholds as low as possible in
the high pile-up environment. Furthermore the level of backgrounds from random
combinations must be maintained or decreased, which might be achieved through
the tracking detectors, the fast timing information in the reconstruction and
the improvement of current selection algorithms based on Machine Learning
tools. Fast timing information to disentangle $pp$-collision points will also
facilitate the analysis of LHCb data and enable the flavour tagging of the
$B^{0}_{s}$ mesons. Further improvements over the projected sensitivities in
this section - especially on $B^{0}\rightarrow\mu^{+}\mu^{-}$ measurements -
might be achieved by improvements on the muon identification and the momentum
resolution, which will have significant impact on the dimuon mass resolution.

###### Acknowledgements.

The authors thank the ATLAS, CMS and LHCb physics working groups for cross
reading the draft. K.-F. Chen is supported by the grant 112-2112-M-002-026 of
National Science and Technology Council, Taiwan. no References

## References

  * Bobeth et al. (2014a) Bobeth, C.; Gorbahn, M.; Hermann, T.; Misiak, M.; Stamou, E.; et al.  $B_{s,d}\rightarrow\ell^{+}\ell^{-}$ in the Standard Model with Reduced Theoretical Uncertainty.  Phys.Rev.Lett. 2014, 112, 101801, [arXiv:hep-ph/1311.0903].  https://doi.org/10.1103/PhysRevLett.112.101801.
  * Bobeth et al. (2014b) Bobeth, C.; Gorbahn, M.; Stamou, E.  Electroweak Corrections to $B_{s,d}\rightarrow\ell^{+}\ell^{-}$.  Phys. Rev. D 2014, 89, 034023, [arXiv:hep-ph/1311.1348].  https://doi.org/10.1103/PhysRevD.89.034023.
  * Hermann et al. (2013) Hermann, T.; Misiak, M.; Steinhauser, M.  Three-loop QCD corrections to $B_{s}\rightarrow\mu^{+}\mu^{-}$.  JHEP 2013, 12, 097, [arXiv:hep-ph/1311.1347].  https://doi.org/10.1007/JHEP12(2013)097.
  * Beneke et al. (2018) Beneke, M.; Bobeth, C.; Szafron, R.  Enhanced electromagnetic correction to the rare B-meson decay $B_{s,d}\rightarrow\mu^{+}\mu^{-}$.  Phys. Rev. Lett. 2018, 120, 011801, [arXiv:hep-ph/1708.09152].  https://doi.org/10.1103/PhysRevLett.120.011801.
  * Beneke et al. (2019) Beneke, M.; Bobeth, C.; Szafron, R.  Power-enhanced leading-logarithmic QED corrections to $B_{q}\rightarrow\mu^{+}\mu^{-}$.  JHEP 2019, 10, 232, [arXiv:hep-ph/1908.07011].  https://doi.org/10.1007/JHEP10(2019)232.
  * Altmannshofer et al. (2012) Altmannshofer, W.; Paradisi, P.; Straub, D.M.  Model-Independent Constraints on New Physics in $b\rightarrow s$ Transitions.  JHEP 2012, 04, 008, [arXiv:hep-ph/1111.1257].  https://doi.org/10.1007/JHEP04(2012)008.
  * Beaujean et al. (2012) Beaujean, F.; Bobeth, C.; van Dyk, D.; Wacker, C.  Bayesian Fit of Exclusive $b\rightarrow s\bar{\ell}\ell$ Decays: The Standard Model Operator Basis.  JHEP 2012, 08, 030, [arXiv:hep-ph/1205.1838].  https://doi.org/10.1007/JHEP08(2012)030.
  * Aoki et al. (2019) Aoki, S.; et al.  FLAG Review 2019, [arXiv:hep-lat/1902.08191].
  * Bazavov et al. (2018) Bazavov, A.; et al.  B- and D-meson leptonic decay constants from four-flavor lattice QCD.  Phys. Rev. D 2018, 98, 074512, [arXiv:hep-lat/1712.09262].  https://doi.org/10.1103/PhysRevD.98.074512.
  * Bussone et al. (2016) Bussone, A.; et al.  Mass of the b quark and B meson decay constants from Nf=2+1+1 twisted-mass lattice QCD.  Phys. Rev. D 2016, 93, 114505, [arXiv:hep-lat/1603.04306].  https://doi.org/10.1103/PhysRevD.93.114505.
  * Dowdall et al. (2013) Dowdall, R.J.; Davies, C.T.H.; Horgan, R.R.; Monahan, C.J.; Shigemitsu, J.  B-meson decay constants from improved lattice nonrelativistic QCD with physical u, d, s, and c quarks.  Phys. Rev. Lett. 2013, 110, 222003, [arXiv:hep-lat/1302.2644].  https://doi.org/10.1103/PhysRevLett.110.222003.
  * Hughes et al. (2018) Hughes, C.; Davies, C.T.H.; Monahan, C.J.  New methods for B meson decay constants and form factors from lattice NRQCD.  Phys. Rev. D 2018, 97, 054509, [arXiv:hep-lat/1711.09981].  https://doi.org/10.1103/PhysRevD.97.054509.
  * De Bruyn et al. (2012) De Bruyn, K.; Fleischer, R.; Knegjens, R.; Koppenburg, P.; Merk, M.; et al.  Probing new physics via the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ effective lifetime.  Phys.Rev.Lett. 2012, 109, 041801, [arXiv:hep-ph/1204.1737].  https://doi.org/10.1103/PhysRevLett.109.041801.
  * Bruyn et al. (2012) Bruyn, K.D.; Fleischer, R.; Knegjens, R.; Koppenburg, P.; Merk, M.; et al.  Branching Ratio Measurements of $B^{0}_{s}$ Decays.  Phys.Rev. 2012, D86, 014027, [arXiv:hep-ph/1204.1735].  https://doi.org/10.1103/PhysRevD.86.014027.
  * Buras (2003) Buras, A.J.  Relations between $\Delta$ M($s$, $d)$ and B($s$, $d)\rightarrow\mu\bar{\mu}$ in models with minimal flavor violation.  Phys. Lett. B 2003, 566, 115–119, [hep-ph/0303060].  https://doi.org/10.1016/S0370-2693(03)00561-6.
  * King et al. (2019) King, D.; Lenz, A.; Rauh, T.  Bs mixing observables and |Vtd/Vts| from sum rules.  JHEP 2019, 05, 034, [arXiv:hep-ph/1904.00940].  https://doi.org/10.1007/JHEP05(2019)034.
  * Buras (2023) Buras, A.J.  Standard Model Predictions for Rare K and B Decays without New Physics Infection.  EPJC 2023, 83, 66, [arXiv:hep-ph/2209.03968].  https://doi.org/10.1140/epjc/s10052-023-11222-6.
  * Amhis et al. (2023) Amhis, Y.; et al.  Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2021.  Phys. Rev. D 2023, 107, 052008, [arXiv:hep-ex/2206.07501].  https://doi.org/10.1103/PhysRevD.107.052008.
  * ATLAS Collaboration (2019) ATLAS Collaboration.  Study of the rare decays of $B^{0}_{s}$ and $B^{0}$ mesons into muon pairs using data collected during 2015 and 2016 with the ATLAS detector.  JHEP 2019, 04, 098, [arXiv:hep-ex/1812.03017].  https://doi.org/10.1007/JHEP04(2019)098.
  * ATLAS Collaboration (2016) ATLAS Collaboration.  Study of the rare decays of $B^{0}_{s}$ and $B^{0}$ into muon pairs from data collected during the LHC Run 1 with the ATLAS detector.  Eur. Phys. J. C 2016, 76, 513, [arXiv:hep-ex/1604.04263].  https://doi.org/10.1140/epjc/s10052-016-4338-8.
  * Hoecker et al. (2007) Hoecker, A.; Speckmayer, P.; Stelzer, J.; Therhaag, J.; von Toerne, E.; Voss, H.  TMVA: Toolkit for Multivariate Data Analysis.  PoS 2007, ACAT, 040, [physics/0703039].
  * Tanabashi et al. (2018) Tanabashi, M.; et al.  Review of Particle Physics.  Phys. Rev. D 2018, 98, 030001.  https://doi.org/10.1103/PhysRevD.98.030001.
  * Amhis et al. (2017) Amhis, Y.; et al.  Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2016.  Eur. Phys. J. C 2017, 77, 895, [arXiv:hep-ex/1612.07233].  https://doi.org/10.1140/epjc/s10052-017-5058-4.
  * (24) Neyman, J.  Outline of a Theory of Statistical Estimation Based on the Classical Theory of Probability.  Phil. Trans. R. Soc. London A, 236 (1937) 333-380.
  * ATLAS Collaboration (2023) ATLAS Collaboration.  Measurement of the $B_{s}^{0}\rightarrow\mu\mu$ Effective Lifetime with the ATLAS Detector.  JHEP 2023, 09, 199, [arXiv:hep-ex/2308.01171].  https://doi.org/10.1007/JHEP09(2023)199.
  * Workman and Others (2022) Workman, R.L.; Others.  Review of Particle Physics.  Progr. Theor. Exp. Phys. 2022, 2022, 083C01 and 2023 update.  https://doi.org/10.1093/ptep/ptac097.
  * Pivk and Le Diberder (2005) Pivk, M.; Le Diberder, F.  sPlot: A statistical tool to unfold data distributions.  Nucl. Instrum. Meth. A 2005, 555, 356–369, [physics/0402083v3].  https://doi.org/https://doi.org/10.1016/j.nima.2005.08.106.
  * CMS Collaboration (2022) CMS Collaboration.  Measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay properties and search for the $B^{0}\rightarrow\mu^{+}\mu^{-}$ decay in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$.  Phys. Lett. B 2022, 842, 137955, [arXiv:hep-ex/2212.10311].  https://doi.org/10.1016/j.physletb.2023.137955.
  * CMS Collaboration (2020) CMS Collaboration.  Measurement of properties of $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decays and search for $B^{0}\rightarrow\mu^{+}\mu^{-}$ with the CMS experiment.  JHEP 2020, 04, 188, [arXiv:hep-ex/1910.12127].  https://doi.org/10.1007/JHEP04(2020)188.
  * Chen and Guestrin (2016) Chen, T.; Guestrin, C.  XGBoost: A Scalable Tree Boosting System 2016.  [arXiv:cs.LG/1603.02754].  https://doi.org/10.1145/2939672.2939785.
  * Pivk and Le Diberder (2005) Pivk, M.; Le Diberder, F.R.  sPlot: a statistical tool to unfold data distributions.  Nucl.Instrum.Meth. 2005, A555, 356–369, [arXiv:physics.data-an/physics/0402083].  https://doi.org/10.1016/j.nima.2005.08.106.
  * Aaij et al. (2021) Aaij, R.; et al.  Precise measurement of the $f_{s}/f_{d}$ ratio of fragmentation fractions and of $B^{0}_{s}$ decay branching fractions.  Phys. Rev. D 2021, 104, 032005, [arXiv:hep-ex/2103.06810].  https://doi.org/10.1103/PhysRevD.104.032005.
  * Read (2002) Read, A.L.  Presentation of search results: The CL(s) technique.  J.Phys. 2002, G28, 2693–2704.  https://doi.org/10.1088/0954-3899/28/10/313.
  * Aaij et al. (2022a) Aaij, R.; et al.  Measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay properties and search for the $B^{0}\rightarrow\mu^{+}\mu^{-}$ and $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}\gamma$ decays.  Phys. Rev. D 2022, 105, 012010, [arXiv:hep-ex/2108.09283].  https://doi.org/10.1103/PhysRevD.105.012010.
  * Aaij et al. (2022b) Aaij, R.; et al.  Analysis of Neutral B-Meson Decays into Two Muons.  Phys. Rev. Lett. 2022, 128, 041801, [arXiv:hep-ex/2108.09284].  https://doi.org/10.1103/PhysRevLett.128.041801.
  * Aaij et al. (2017) Aaij, R.; et al.  Measurement of the $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\rightarrow\mu^{+}\mu^{-}$ decays.  Phys. Rev. Lett. 2017, 118, 191801, [arXiv:hep-ex/1703.05747].  https://doi.org/10.1103/PhysRevLett.118.191801.
  * Archilli et al. (2013) Archilli, F.; et al.  Performance of the muon identification at LHCb.  JINST 2013, 8, P10020, [arXiv:physics.ins-det/1306.0249].  https://doi.org/10.1088/1748-0221/8/10/P10020.
  * Zyla et al. (2020) Zyla, P.A.; et al.  Review of Particle Physics.  PTEP 2020, 2020, 083C01.  https://doi.org/10.1093/ptep/ptaa104.
  * Skwarnicki (1986) Skwarnicki, T.  A study of the radiative cascade transitions between the Upsilon-prime and Upsilon resonances.  PhD thesis, Institute of Nuclear Physics, Krakow, 1986.  DESY-F31-86-02.
  * ATL (2020) Combination of the ATLAS, CMS and LHCb results on the $B^{0}_{(s)}\rightarrow\mu^{+}\mu^{-}$ decays.  Technical report, CERN, 2020.  ATLAS-CONF-2020-049, CMS-PAS-BPH-20-003, LHCb-CONF-2020-002.
  * Sirunyan et al. (2019) Sirunyan, A.M.; et al.  Measurement of properties of B${}^{0}_{\mathrm{s}}\rightarrow\mu^{+}\mu^{-}$ decays and search for B${}^{0}\rightarrow\mu^{+}\mu^{-}$ with the CMS experiment 2019.  [arXiv:hep-ex/1910.12127].
  * Khachatryan et al. (2015) Khachatryan, V.; et al.  Observation of the rare $B^{0}_{s}\rightarrow\mu^{+}\mu^{-}$ decay from the combined analysis of CMS and LHCb data.  Nature 2015, 522, 68–72, [arXiv:hep-ex/1411.4413].  https://doi.org/10.1038/nature14474.
  * ATL (2022) Snowmass White Paper Contribution: Physics with the Phase-2 ATLAS and CMS Detectors 2022.
  * LHC (2018) Physics case for an LHCb Upgrade II: Opportunities in flavour physics, and beyond, in the HL-LHC era 2018.
  * Buras et al. (2013) Buras, A.J.; Fleischer, R.; Girrbach, J.; Knegjens, R.  Probing New Physics with the $B_{s}\rightarrow{\mu}+{\mu}-$ Time-Dependent Rate.  JHEP 2013, 07, 077, [arXiv:hep-ph/1303.3820].  https://doi.org/10.1007/JHEP07(2013)077.