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

TL;DR

This paper reports on the measurement of rare $B_s^0$ and $B^0$ meson decays into muon pairs using ATLAS data from 2015-2016, providing results consistent with the Standard Model within uncertainties.

Contribution

First measurement combining 2015-2016 ATLAS data for these rare decays, improving precision and confirming Standard Model predictions.

Findings

Measured ${ m B}(B_s^0 o \mu^+\mu^-)$ = (3.2^{+1.1}_{-1.0}) x 10^{-9}

Set upper limit ${ m B}(B^0 o \mu^+\mu^-) < 4.3 imes 10^{-10}$ at 95% CL

Combined results with Run 1 data are consistent with Standard Model

Abstract

A study of the decays $B^0_s \to \mu^+\mu^-$ and $B^0 \to \mu^+\mu^-$ has been performed using 26.3 fb$^{-1}$ of 13 TeV LHC proton-proton collision data collected with the ATLAS detector in 2015 and 2016. Since the detector resolution in $\mu^+\mu^-$ invariant mass is comparable to the $B_{s}^{0}$-$B^0$ mass difference, a single fit determines the signal yields for both decay modes. This results in a measurement of the branching fraction ${\cal B}(B^0_s \to \mu^+\mu^-) = \left( 3.2^{+1.1}_{-1.0} \right) \times 10^{-9}$ and an upper limit ${\cal B}(B^0 \to \mu^+\mu^-) < 4.3 \times 10^{-10}$ at 95% confidence level. The result is combined with the Run 1 ATLAS result, yielding ${\cal B}(B^0_s \to \mu^+\mu^-) = \left( 2.8^{+0.8}_{-0.7} \right) \times 10^{-9}$ and ${\cal B}(B^0 \to \mu^+\mu^-) < 2.1 \times 10^{-10}$ at 95% confidence level. The combined result is consistent with the Standard Model prediction within 2.4 standard deviations in the ${\cal B}(B^0 \to \mu^+\mu^-)$-${\cal B}(B^0_s \to \mu^+\mu^-)$ plane.

Tables

Table 1. Table 1: Description of the 15 15 15 input variables used in a BDT classifier to discriminate between signal and continuum background. When the BDT classifier is applied to B + → J / ψ ​ K + → superscript 𝐵 𝐽 𝜓 superscript 𝐾 B^{+}\to J/\psi\,K^{+}  and B s 0 → J / ψ ​ ϕ → subscript superscript 𝐵 0 𝑠 𝐽 𝜓 italic-ϕ B^{0}_{s}\to J/\psi\,\phi  candidates, the variables related to the decay products of the B 𝐵 B mesons refer only to the muons from the decay of the J / ψ 𝐽 𝜓 J/\psi . Horizontal lines separate the classifications into groups (a), (b) and (c) respectively, as described in the text. For category (c), additional tracks are required to have p T subscript 𝑝 T p_{\text{T}} >500 MeV .

Variable	Description
$p_{\text{T}}^B$	Magnitude of the $B$ candidate transverse momentum ${\overrightarrow{p_{\text{T}}}}^B$.
${\chi }_{\mathrm{PV},\mathrm{DV}xy}^2$	Compatibility of the separation $\overrightarrow{\Delta x}$ between production (i.e. associated PV) and decay (DV) vertices in the transverse projection: ${\overrightarrow{\Delta x}}_{\text{T}}\cdot {\Sigma }_{{\overrightarrow{\Delta x}}_{\text{T}}}^{-1}\cdot {\overrightarrow{\Delta x}}_{\text{T}}$, where ${\Sigma }_{{\overrightarrow{\Delta x}}_{\text{T}}}$ is the covariance matrix.
$\Delta R_{\mathrm{flight}}$	Three-dimensional angular distance between ${\overrightarrow{p}}^B$ and $\overrightarrow{\Delta x}$: $\sqrt{{\alpha }_{2D}{}_{}{}^2+{\left(\Delta \eta \right)}^2}$
$\left|{\alpha }_{2D}\right|$	Absolute value of the angle in the transverse plane between ${\overrightarrow{p_{\text{T}}}}^B$ and ${\overrightarrow{\Delta x}}_{\text{T}}$.
$L_{xy}$	Projection of ${\overrightarrow{\Delta x}}_{\text{T}}$ along the direction of ${\overrightarrow{p}}_{\text{T}}^B$: $\left({\overrightarrow{\Delta x}}_{\text{T}}\cdot {\overrightarrow{p_{\text{T}}}}^B\right)/\left|{\overrightarrow{p_{\text{T}}}}^B\right|$.
IP${}_{}{}^{3D}{}_B{}^{}$	Three-dimensional impact parameter of the $B$ candidate to the associated PV.
DOCAμμ	Distance of closest approach (DOCA) of the two tracks forming the $B$ candidate (three-dimensional).
$\Delta {\varphi }_{\mu \mu }$	Azimuthal angle between the momenta of the two tracks forming the $B$ candidate.
${\left|d_0\right|}^{\max }$-sig.	Significance of the larger absolute value of the impact parameters to the PV of the tracks forming the $B$ candidate, in the transverse plane.
${\left|d_0\right|}^{\min }$-sig.	Significance of the smaller absolute value of the impact parameters to the PV of the tracks forming the $B$ candidate, in the transverse plane.
$P_{\text{L}}^{\min }$	The smaller of the projected values of the muon momenta along ${\overrightarrow{p_{\text{T}}}}^B$.
$I_{0.7}$	Isolation variable defined as ratio of $\left|{\overrightarrow{p_{\text{T}}}}^B\right|$ to the sum of $\left|{\overrightarrow{p_{\text{T}}}}^B\right|$ and the transverse momenta of all additional tracks contained within a cone of size $\Delta R=\sqrt{{\left(\Delta \varphi \right)}^2+{\left(\Delta \eta \right)}^2}=0.7$ around the $B$ direction. Only tracks matched to the same PV as the $B$ candidate are included in the sum.
DOCAxtrk	DOCA of the closest additional track to the decay vertex of the $B$ candidate. Only tracks matched to the same PV as the $B$ candidate are considered.
$N_{\mathrm{xtrk}}^{\mathrm{close}}$	Number of additional tracks compatible with the decay vertex (DV) of the $B$ candidate with . Only tracks matched to the same PV as the $B$ candidate are considered.
${\chi }_{\mu ,\mathrm{xPV}}^2$	Minimum ${\chi }^2$ for the compatibility of a muon in the $B$ candidate with any PV reconstructed in the event.

Table 2. Table 2: Summary of the uncertainties in R ε subscript 𝑅 𝜀 R_{\varepsilon} .

Source	Contribution [%]
Statistical	0.8
Kinematic reweighting (DDW)	0.8
Muon trigger and reconstruction	1.0
BDT input variables	3.2
Kaon tracking efficiency	1.5
Pile-up reweighting	0.6

Table 3. Table 3: Breakdown of the expected systematic uncertainties in ℬ ℬ \cal B ( B ( s ) 0 → μ + ​ μ − → subscript superscript 𝐵 0 𝑠 superscript 𝜇 superscript 𝜇 B^{0}_{(s)}\to\mu^{+}\mu^{-} ). The measurements are dominated by statistical uncertainty, followed by the systematic uncertainty from the fit. The latter is dominated by contributions from the mass scale uncertainty and the parameterisation of the b → μ + ​ μ − ​ X → 𝑏 superscript 𝜇 superscript 𝜇 𝑋 b\to\mu^{+}\mu^{-}X background. The statistical uncertainties reported here are obtained from the maximisation of the fit likelihood and are meant only as a reference for the relative scale uncertainties.

Source	$B_s^0$ [%]	$B^0$ [%]
$f_s/f_d$	5.1	-
$B^{+}$ yield	4.8	4.8
$R_{\epsilon }$	4.1	4.1
$ℬ\left(B^{+}\to J/\psi K^{+}\right)\times ℬ\left(J/\psi \to {\mu }^{+}{\mu }^{-}\right)$	2.9	2.9
Fit systematic uncertainties	8.7	65
Stat. uncertainty (from likelihood est.)	27	150

Equations

$$\mathcal{B}(B_{(s)}^0\to {\mu }^{+}{\mu }^{-})$$

$$\mathcal{B}(B_s^0\to {\mu }^{+}{\mu }^{-})=\left(3.21_{-0.91-0.30}^{+0.96+0.49}\right)\times 10^{-9}=\left(3.2_{-1.0}^{+1.1}\right)\times 10^{-9}.$$

$$\mathcal{B}(B^0\to {\mu }^{+}{\mu }^{-})<4.3\text{leavevmode}\times \text{leavevmode}10^{-10}.$$

$$\mathcal{B}(B_s^0\to {\mu }^{+}{\mu }^{-})=\left(2.8\pm 0.7\right)\times 10^{-9},$$

$$\mathcal{B}(B^0\to {\mu }^{+}{\mu }^{-})=\left(-1.9\pm 1.6\right)\times 10^{-10}.$$

$$\mathcal{B}(B_s^0\to {\mu }^{+}{\mu }^{-})=\left(2.8_{-0.7}^{+0.8}\right)\times 10^{-9}.$$

$$\mathcal{B}(B^0\to {\mu }^{+}{\mu }^{-})<2.1\text{leavevmode}\times \text{leavevmode}10^{-10}.$$

Full text

\PreprintIdNumber

CERN-EP-2018-291

\AtlasTitleStudy 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 \AtlasVersion2.9 \AtlasAbstract A study of the decays $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$ has been performed using $26.3$ fb*-1* of $13$ TeV LHC proton–proton collision data collected with the ATLAS detector in 2015 and 2016. Since the detector resolution in $\mu^{+}\mu^{-}$ invariant mass is comparable to the $B_{s}^{0}$–$B^{0}$ mass difference, a single fit determines the signal yields for both decay modes. This results in a measurement of the branching fraction ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(3.2^{+1.1}_{-1.0}\right)\times 10^{-9}$ and an upper limit ${\cal B}(B^{0}\to\mu^{+}\mu^{-})<4.3\times 10^{-10}$ at $95\%$ confidence level. The result is combined with the Run 1 ATLAS result, yielding ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(2.8^{+0.8}_{-0.7}\right)\times 10^{-9}$ and ${\cal B}(B^{0}\to\mu^{+}\mu^{-})<2.1\times 10^{-10}$ at 95% confidence level. The combined result is consistent with the Standard Model prediction within 2.4 standard deviations in the ${\cal B}(B^{0}\to\mu^{+}\mu^{-})$-${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})$ plane.

\AtlasRefCodeBPHY-2018-09 \AtlasJournalJournal of High Energy Physics \AtlasJournalRefJHEP 04 (2019) 098 \AtlasDOI10.1007/JHEP04(2019)098 \AtlasCoverSupportingNoteANA-BPHY-2018-09-INT1https://cds.cern.ch/record/2311593 \AtlasCoverCommentsDeadline31 Oct 2018 \AtlasCoverAnalysisTeamA. Cerri, I. Ibragimov, A. Grummer, F. Tresoldi, U. de Sanctis, A. Campoverde, W. Walkowiak, P. Buchholz, S. Seidel, S.Yu. Sivoklokov \AtlasCoverEdBoardMemberJulie Kirk, Bill Murray (chair), Sandro Palestini \AtlasCoverEgroupEditorsatlas-BPHY-2018-09-editors@cern.ch \AtlasCoverEgroupEdBoardatlas-BPHY-2018-09-editorial-board@cern.ch

1 Introduction

Flavour-changing neutral-current processes are highly suppressed in the Standard Model (SM). The branching fractions of the decays $B^{0}_{(s)}\to\mu^{+}\mu^{-}$  are, in addition, helicity suppressed in the SM, and are predicted to be ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=(3.65\pm 0.23)\times 10^{-9}$ and ${\cal B}(B^{0}\to\mu^{+}\mu^{-})=(1.06\pm 0.09)\times 10^{-10}$ [1]. The smallness and precision of these predicted branching fractions provide a favourable environment for observing contributions from new physics: significant deviations from SM predictions could arise in models involving non-SM heavy particles, such as those predicted in the Minimal Supersymmetric Standard Model [2, 3, 4, 5, 6] and in extensions such as Minimal Flavour Violation [7, 8], Two-Higgs-Doublet Models [6], and others [9, 10]. The CMS and LHCb collaborations reported the observation of $B^{0}_{s}\to\mu^{+}\mu^{-}$ [11, 12] and evidence of $B^{0}\to\mu^{+}\mu^{-}$ with combined values ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(2.8^{+0.7}_{-0.6}\right)\times 10^{-9}$ and ${\cal B}(B^{0}\to\mu^{+}\mu^{-})=\left(3.9^{+1.6}_{-1.4}\right)\times 10^{-10}$ [13]. The LHCb Collaboration updated its Run 1 result with part of the data collected in Run 2, measuring ${\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) and determining ${\cal B}(B^{0}\to\mu^{+}\mu^{-})<3.4\times 10^{-10}$ at $95\%$ confidence level (CL) [14]. ATLAS has measured, with the Run 1 dataset [15], ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(0.9^{+1.1}_{-0.8}\right)\times 10^{-9}$, while setting an upper limit on $B^{0}\to\mu^{+}\mu^{-}$ of $4.2\times 10^{-10}$ at $95\%$ CL.

This paper reports the result of a search for $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$  decays performed using $pp$ collision data corresponding to an effective integrated luminosity of 26.3 fb*-1*, collected at $13$ TeV centre-of-mass energy during the first two years of the LHC Run 2 data-taking period using the ATLAS detector. The analysis strategy follows the approach employed in the previous ATLAS measurement [15], but uses data collected with one rather than three separate sets of trigger thresholds, applies standard ATLAS muon selection criteria (Section 3) rather than a dedicated muon multivariate discriminant, and employs an improved statistical treatment of the result (Section 10).

The notation used throughout the paper refers to the combination of processes and their charge-conjugates, unless otherwise specified. The $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$ branching fractions are measured relative to the reference decay mode $B^{+}\to J/\psi(\to\mu^{+}\mu^{-})K^{+}$ which is abundant and has a well-measured branching fraction ${\cal B}(B^{+}\to J/\psi\,K^{+})\times{\cal B}(J/\psi\to\mu^{+}\mu^{-})$. The $B^{0}\to\mu^{+}\mu^{-}$ ($B^{0}_{s}\to\mu^{+}\mu^{-}$) branching fraction can be extracted as:

[TABLE]

where $N_{d}$ ($N_{s}$) is the $B^{0}\to\mu^{+}\mu^{-}$ ($B^{0}_{s}\to\mu^{+}\mu^{-}$) signal yield, $N_{J/\psi K^{+}}$ is the $B^{+}\to J/\psi\,K^{+}$ reference channel yield, $\varepsilon_{\mu^{+}\mu^{-}}$ and $\varepsilon_{J/\psi K^{+}}$ are the corresponding values of acceptance times efficiency (measured in fiducial regions defined in Section 9), and $f_{u}/f_{d}$ ($f_{u}/f_{s}$) is the ratio of the hadronisation probabilities of a $b$-quark into $B^{+}$ and $B^{0}$ ($B_{s}^{0}$). In the quantity $\mathcal{D}_{\mathrm{ref}}=N_{J/\psi K^{+}}\times(\varepsilon_{\mu^{+}\mu^{-}}/\varepsilon_{J/\psi K^{+}})$, the $\varepsilon$ ratio takes into account relative differences in efficiencies, integrated luminosities and the trigger selections used for the signal and the reference modes. Signal and reference channel events are selected with similar dimuon triggers. One half of the reference channel sample is used to determine the normalisation and the other half is used to tune the kinematic distributions of simulated events.

The event selection uses variables related to the $B$ candidate decay time, thus introducing a dependence of the efficiency on the signal lifetime. The relation between the measured branching fraction and the corresponding value at production is established assuming the decay time distribution predicted in the SM, where the decay occurs mainly through the heavy eigenstate $B^{0}_{(s),\mathrm{H}}$ of the $B^{0}_{(s)}$–$\overline{B^{0}_{(s)}}$ system. Some models of new physics [16, 17] predict modifications to the decay time distribution of $B^{0}_{s}\to\mu^{+}\mu^{-}$ and a comparison with the experimental result requires a correction to the ratio of the time-integrated efficiencies entering $\mathcal{D}_{\mathrm{ref}}$.

The ATLAS inner tracking system, muon spectrometer and, for efficient identification of muons, also the calorimeters, are used to reconstruct and select the event candidates. Details of the detector, trigger, data sets, and preliminary selection criteria are discussed in Sections 2 and 3. A blind analysis was performed in which data in the dimuon invariant mass range from $5166$ to $5526$ MeV were removed until the procedures for event selection and the details of signal yield extraction were completely defined. Section 4 introduces the three main categories of background. Section 5 describes the strategy used to reduce the probability of hadron misidentification. The final sample of candidates is selected using a multivariate classifier, designed to enhance the signal relative to the dominant dimuon background component, as discussed in Section 6. Checks on the distributions of the variables used in the multivariate classifier are summarised in Section 7. They are based on the comparison of data and simulation for dimuon events, for $B^{+}\to J/\psi\,K^{+}$ candidates and for events selected as $B^{0}_{s}\to J/\psi\,\phi\to\mu^{+}\mu^{-}K^{+}K^{-}$, which provide an additional validation of the procedures used in the analysis. Section 8 details the fit procedure used to extract the yield of $B^{+}\to J/\psi\,K^{+}$ events. The determination of the ratio of efficiencies in the signal and the reference channels is presented in Section 9. Section 10 describes the extraction of the signal yield, obtained with an unbinned maximum-likelihood fit performed on the dimuon invariant mass distribution. In this fit, events are separated into classifier intervals to maximise the fit sensitivity. The results for the branching fractions $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) and $\cal B$($B^{0}\to\mu^{+}\mu^{-}$) are reported in Section 11 and combined with the full Run 1 results in Section 12.

2 ATLAS detector, data and simulation samples

The ATLAS detector111ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point. The $z$-axis is along the beam pipe, the $x$-axis points to the centre of the LHC ring and the $y$-axis points upward. Cylindrical coordinates $(r,\phi)$ are used in the transverse plane, $r$ being the distance from the origin and $\phi$ being the azimuthal angle around the beam pipe. The pseudorapidity $\eta$ is defined as $\eta=-\ln[\tan(\theta/2)]$ where $\theta$ is the polar angle. consists of three main components: an inner detector (ID) tracking system immersed in a $2$ T axial magnetic field, surrounded by electromagnetic and hadronic calorimeters and by the muon spectrometer (MS). A full description can be found in Ref. [18], complemented by Ref. [19] for details about the new innermost silicon pixel layer that was installed for Run 2.

This analysis is based on the Run 2 data recorded in 2015 and 2016 from $pp$ collisions at the LHC at $\sqrt{s}=13$ TeV. Data used in the analysis were recorded during stable LHC beam periods. Data quality requirements were imposed, notably on the performance of the MS, ID and calorimeter systems. The total integrated luminosity collected by ATLAS in this period is $36.2$ fb*-1* with an uncertainty of 2.1%. These values are determined using a methodology similar to that detailed in Ref. [20], based on calibration of the luminosity scale using $x$–$y$ beam-separation scans, and use the LUCID-2 detector [21] for the baseline luminosity measurement. The total effective integrated luminosity used in this analysis - accounting for trigger prescales - amounts to $26.3$ fb*-1* for the signal and $15.1$ fb*-1* for the reference channel.

Samples of simulated Monte Carlo (MC) events are used for training and validation of the multivariate analyses, for the determination of the efficiency ratios, and for developing the procedure used to determine the signal. Exclusive MC samples were produced for the signal channels $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$, the reference channel $B^{+}\to J/\psi\,K^{+}$ ($J/\psi\to\mu^{+}\mu^{-}$), and the control channel $B^{0}_{s}\to J/\psi\,\phi$ ($J/\psi\to\mu^{+}\mu^{-}$, $\phi\to K^{+}K^{-}$). In addition, background studies employ MC samples of inclusive $b\bar{b}\to\mu^{+}\mu^{-}X$ decays, exclusive samples of $B_{s}^{0}\to K^{-}\mu^{+}\nu$, $B^{0}\to\pi^{-}\mu^{+}\nu$, $\Lambda_{b}\to p\mu^{-}\overline{\nu}$, $B^{0}_{(s)}\to hh^{\prime}$ decays with $h^{(\prime)}$ being a charged pion or kaon, and inclusive decays $B\to J/\psi X$ as well as the exclusive $B^{+}\to J/\psi\,\pi^{+}$ decay.

Most of the dimuon candidates in the data sample originate from the decays of hadrons produced in the hadronisation of $b\bar{b}$ pairs. The inclusive $b\bar{b}\to\mu^{+}\mu^{-}X$ MC sample used to describe this background requires the presence of two muons in the final state, with both muons originating from the $b\bar{b}$ decay chain. The size of this sample is equivalent to roughly three times the integrated luminosity of the data.

The MC samples were generated with Pythia 8 [22]. The ATLAS detector and its response were simulated using Geant4 [23, 24]. Additional $pp$ interactions in the same and nearby bunch crossings (pile-up) are included in the simulation. Muon reconstruction and triggering efficiencies are corrected in the simulated samples using data-driven scale factors. The scale factors for the trigger efficiencies are obtained by comparing data and simulation efficiencies determined with a tag-and-probe method on inclusive prompt and non-prompt $J/\psi$ candidates. This procedure yields scale factors as a function of the muon transverse momentum and pseudorapidity, which are applied throughout the analysis [25]. Reconstruction and selection efficiencies are obtained from simulation and similarly corrected according to data-driven comparisons. In addition to these efficiency corrections, simulated events are reweighted to reproduce the pile-up multiplicity observed in data, and according to the equivalent integrated luminosity associated with each trigger selection.

Using the iterative reweighting method described in Ref. [26], the simulated samples of the exclusive decays considered are adjusted with two-dimensional data-driven weights (DDW) to correct for the differences between simulation and data observed in the $B$ meson transverse momentum and pseudorapidity distributions. DDW obtained from $B^{+}\to J/\psi\,K^{+}$ decays are used to correct the simulation samples in the signal and reference channels. DDW obtained from the $B^{0}_{s}\to J/\psi\,\phi$ control channel are found to agree with those from $B^{+}\to J/\psi\,K^{+}$, showing that the same corrections are applicable to $B_{s}^{0}$ and $B^{0}$ decays.

Residual differences between data and simulation studied in the $B^{+}\to J/\psi\,K^{+}$ and $B^{0}_{s}\to J/\psi\,\phi$ signals are treated as sources of systematic uncertainty in the evaluation of the signal efficiency, as discussed in Section 9. The only exception to this treatment is the $B$ meson isolation ($I_{0.7}$ in Section 6 and Table 1), where residual differences are used to reweight the signal MC events and the corresponding uncertainties are propagated to account for residual systematic uncertainty effects.

Similarly to the exclusive decays, the kinematic distributions of the inclusive $b\bar{b}\to\mu^{+}\mu^{-}X$ MC sample are reweighted with corrections obtained from the dimuon invariant mass sidebands in data.

3 Data selection

For data collected during the LHC Run 2, the ATLAS detector uses a two-level trigger system, consisting of a hardware-based first-level trigger and a software-based high-level trigger. A first-level dimuon trigger [27] selects events requiring that one muon has $p_{\text{T}}>4$ GeV and the other has $p_{\text{T}}>6$ GeV. A full track reconstruction of the muon candidates is performed by the high-level trigger, where an additional loose selection is imposed on the dimuon invariant mass $m_{\mu\mu}$, accepting candidates in the range 4 GeV to 8.5 GeV. Due to the increased pile-up in 2016 data, an additional selection was added at this trigger stage, requiring the vector from the primary vertex to the dimuon vertex to have a positive component ($L_{xy}$) along the dimuon’s transverse momentum direction. The effect of this selection is accounted for in the analysis but has no consequence since stricter requirements are applied in the full event selection (see Section 6).

The signal channel, the reference channel $B^{+}\to J/\psi\,K^{+}$ and the control channel $B_{s}^{0}\to J/\psi\phi$ were selected with trigger prescale factors that vary during the data-taking period. In the 36.2 fb*-1* of data analysed, the prescaling of the trigger approximately averages to a reduction by a factor 1.4, giving an effective integrated luminosity for the signal sample of 26.3 fb*-1*, while for the reference and control channels 15.1 fb*-1* were collected due to an effective prescale of 2.4. These effects are taken into account in the extraction of the signal branching fraction, through the $\varepsilon$ factors in Eq. (1).

Using information from the full offline reconstruction, a preliminary selection is performed on candidates for $B^{0}_{(s)}\to\mu^{+}\mu^{-}$, $B^{+}\to J/\psi\,K^{+}\to\mu^{+}\mu^{-}K^{+}$ and $B_{s}^{0}$ $\to J/\psi\phi\to\mu^{+}\mu^{-}K^{+}K^{-}$ decays. In the ID system, muon candidates are required to have at least one measured hit in the pixel detector and two measured hits in the semiconductor tracker. They are also required to be reconstructed in the MS, and to have $|\eta|<2.5$. The offline muon pair must pass the $p_{\text{T}}>4$ GeV and $p_{\text{T}}>6$ GeV requirements imposed by the trigger. Furthermore, the muon candidates are required to fulfil tight muon quality criteria [28]; this requirement is relaxed to loose for the hadron misidentification studies in Section 5. Kaon candidates must satisfy similar requirements in the ID, except for a looser requirement of $p_{\text{T}}>1$ GeV.

The computed $B$ meson properties are based on a decay vertex fitted to two, three or four tracks, depending on the decay process to be reconstructed. The $B$ candidates are required to have a $\chi^{2}$ per degree of freedom below 6 for the fit to the $B$ vertex, and below 10 for the fit to the $J/\psi\to\mu^{+}\mu^{-}$ vertex. The selections $2915<m(\mu^{+}\mu^{-})<3275$ MeV and $1005<m(K^{+}K^{-})<1035$ MeV are applied to the $J/\psi\to\mu^{+}\mu^{-}$ and the $\phi\to K^{+}K^{-}$ vertices, respectively. In the fits to the $B^{+}\to J/\psi\,K^{+}$ and $B^{0}_{s}\to J/\psi\,\phi$ channels, the reconstructed dimuon mass is constrained to the world average $J/\psi$ mass [29].

Reconstructed $B$ candidates are retained if they satisfy $p_{\text{T}}^{B}>8.0$ GeV and $|\eta^{B}|<2.5$. The invariant mass of each $B$ candidate is calculated using muon trajectories measured by combining the information from the ID and MS to improve upon the mass resolution obtained from ID information only [30].

The invariant mass range considered for the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ decay starts at $4766$ MeV and is $1200$ MeV wide. Within this range a $360$-MeV-wide signal region is defined, starting at $5166$ MeV. The remainder of the range defines the upper and lower mass sidebands of the analysis.

For the reference and control channels, the mass range considered is $4930$–$5630$ ($5050$–$5650$) MeV for $B^{+}\to J/\psi\,K^{+}$ ($B^{0}_{s}\to J/\psi\,\phi$), where $5180$–$5380$ ($5297$–$5437$) MeV is the peak region and higher and lower mass ranges comprise the mass sidebands used for background subtraction.

The coordinates of primary vertices (PV) are obtained from charged-particle tracks not used in the decay vertices, and that are constrained to the luminous region of the colliding beams in the transverse plane. The matching of a $B$ candidate to a PV is made by extrapolating the candidate trajectory to the point of closest approach to the beam axis, and choosing the PV with the smallest distance along $z$. Simulation shows that this method matches the correct vertex with a probability above $99\%$ for all relevant pile-up conditions.

To reduce the large background in the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ channel before applying the final selection based on multivariate classifiers, a loose collinearity requirement is applied between the momentum of the $B$ candidate ($\overrightarrow{p}^{B}$) and the vector from the PV to the decay vertex ($\overrightarrow{\Delta x}$). The absolute value of the azimuthal angle $\alpha_{2\mathrm{D}}$ between these two vectors is required to be smaller than 1.0 radians. The combination $\Delta R_{\mathrm{flight}}=\sqrt{\alpha_{2\mathrm{D}}{}^{2}+(\Delta\eta)^{2}}$, where $\Delta\eta$ is the difference in pseudorapidity, is required to satisfy $\Delta R_{\mathrm{flight}}<1.5$.

After the preliminary selection, approximately $3.5\times 10^{6}$ candidates are found in the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ fit region, with about $1.0\times 10^{6}$ falling in the blinded range [5166, 5526] MeV.

4 Background composition

The background to the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ signal originates from three main sources:

Continuum background, the dominant combinatorial component, which consists of muons originating from uncorrelated hadron decays and is characterised by a weak dependence on the dimuon invariant mass;

Partially reconstructed decays, where one or more of the final-state particles ($X$) in a $b$-hadron decay is not reconstructed, causing these candidates to accumulate in the low dimuon invariant mass sideband (this background includes a significant contribution from semileptonic decays where one of the muons is a misidentified hadron, discussed below);

Peaking background, due to $B^{0}_{(s)}\to hh^{\prime}$ decays, with both hadrons misidentified as muons.

The continuum background consists mainly of muons produced independently in the fragmentation and decay chains of a $b$-quark and a $\bar{b}$-quark. It is studied in the signal mass sidebands, and it is found to be well described by the inclusive $b\bar{b}\to\mu^{+}\mu^{-}X$ MC sample.

The partially reconstructed decays consist of several topologies: (a) same-side combinatorial background from decay cascades ($b\to c\mu^{-}\overline{\nu}\to s(d)\mu^{+}\mu^{-}\nu\overline{\nu}$); (b) same-vertex background from $B$ decays containing a muon pair (e.g. $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ or $B\to J/\psi\mu^{+}X\to\mu^{+}\mu^{-}\mu^{+}X$); (c) $B_{c}^{+}$ decays (e.g. $B_{c}^{+}\to J/\psi\mu^{+}\nu\to\mu^{+}\mu^{-}\mu^{+}\nu$); (d) semileptonic $b$-hadron decays where a final-state hadron is misidentified as a muon. The remainder of this paper implicitly excludes categories (c) and (d) when referring to partially reconstructed or $b\to\mu^{+}\mu^{-}X$ decays, since these categories are treated separately.

The $b\bar{b}\to\mu^{+}\mu^{-}X$ MC sample is used to investigate the background composition after the analysis selection. All backgrounds in this sample have a dimuon invariant mass distribution mainly below the mass range considered in this analysis, with a high-mass tail extending through the signal region. The simulation does not contemplate sources other than muons from $b\bar{b}$ decays: $c\bar{c}$ and prompt contributions are not included222These sources are suppressed by the final analysis selections introduced in Section 6. Potential residual contributions are found to be consistent with the continuum background models used in the final fit (Section 10.1), where systematic uncertainties for background model inconsistencies between data and MC are taken into account (Section 10.3).. All possible origins of two muons in the $b\bar{b}$ decay tree are, however, analysed, after classification into the mutually exclusive continuum and partially reconstructed categories described above. This sample is used only to identify suitable functional models for the corresponding background components, and as a benchmark for these models. No shape or normalisation constraints are derived from this simulation. This makes the analysis largely insensitive to mismatches between background simulation and data.

The semileptonic decays with final-state hadrons misidentified as muons consist mainly of three-body charmless decays $B^{0}\to\pi^{-}\mu^{+}\nu$, $B^{0}_{s}\to K^{-}\mu^{+}\nu$ and $\Lambda_{b}\to p\mu^{-}\overline{\nu}$ in which the tail of the invariant mass distribution extends into the signal region. Due to branching fractions of the order of $10^{-6}$, this background is not large, and is further reduced by the muon identification requirements, discussed in Section 5. The MC invariant mass distributions of these partially reconstructed decay topologies are shown together with the SM signal predictions in Figure 1(a) after applying the preliminary selection criteria described in Section 3.

Finally, the peaking background is due to $B^{0}_{(s)}\to hh^{\prime}$ decays containing two hadrons misidentified as muons. The distributions in Figure 1(b), obtained from simulation, show that these decays populate the signal region. This component is further discussed in Section 5.

5 Hadron misidentification

In the preliminary selection, muon candidates are formed from the combination of tracks reconstructed independently in the ID and MS. The performance of the muon reconstruction in ATLAS is presented in Ref. [28]. Additional studies were performed to evaluate the amount of background related to hadrons erroneously identified as muons.

Detailed simulation studies were performed for the $B^{0}_{(s)}\to hh^{\prime}$ channel with a full Geant4-based simulation [23] of all systems of the ATLAS detector. The vast majority of background events from particle misidentification are due to decays in flight of kaons and pions, in which the muon receives most of the energy of the parent meson. Hence this background is generally related to true muons measured in the MS, but not produced promptly in the decay of a $B$ meson.

The muon candidate is required to pass tight muon requirements in the preliminary selection, which are based on the profile of energy deposits in the calorimeters as well as on tighter ID–MS matching criteria than those used for the loose requirements. Two-body $B$ decays in control regions show that tight selections have, relative to the loose counterpart, an average hadron misidentification probability reduced by a factor 0.39 with a muon reconstruction efficiency of 90%. The resulting final value of the misidentification probability is 0.08% for kaons and 0.1% for pions. Efficiencies and fake rates are relative to the analysis preselections, including tracking but excluding any muon requirement.

The background due to $B^{0}_{(s)}\to hh^{\prime}$, with double misidentification of $hh^{\prime}$ as $\mu^{+}\mu^{-}$, has a reconstructed invariant mass distribution that peaks at 5240 MeV, close to the $B^{0}$ mass, and is effectively indistinguishable from the $B^{0}$ signal (see Figure 1(b)). The expected number of peaking-background events can be estimated in a way analogous to that for the signal, from the number of observed $B^{+}\to J/\psi\,K^{+}$ events using Eq. (1), after taking into account the expected differences from muon identification variables and trigger selections. World average [29] values for the branching fractions of $B^{0}$ and $B_{s}^{0}$ into $K\pi$, $KK$ and $\pi\pi$ are used, together with the hadron misidentification probabilities obtained from simulation. This results in $2.7\pm 1.3$ total expected peaking-background events, after the reference multivariate selection.333This selection, corresponding to 54% signal efficiency, was also applied to derive all other quantities quoted in this Section and includes a selection against the $\mu^{+}\mu^{-}$ continuum background based on the BDT discussed in Section 6.

When selecting loose muons and inverting the additional requirements imposed in the tight muon selection, the number of events containing real muons is substantially reduced, while the number of peaking-background events is approximately two times larger than in the sample obtained with the nominal selection. A fit to data for this background-enhanced sample returns $6.8\pm 3.7$ events, which translates into a peaking-background yield in the signal region of $2.9\pm 2.0$ events when taking into account the relative rejection of the muon quality selections. The predicted yield is in good agreement with the simulation.

Besides the peaking background, the tight muon selection also reduces the semileptonic contributions with a single misidentified hadron. Simulation yields $30\pm 3$ events expected from $B^{0}\to\pi^{-}\mu^{+}\nu$ and $B_{s}^{0}\to K^{-}\mu^{+}\nu$ in the final sample, with a distribution kinematically constrained to be mostly below the signal region. The $\Lambda_{b}\to p\mu^{-}\overline{\nu}$ contribution is negligible due to the smaller production cross section and the low rate at which protons fake muons.

6 Continuum background reduction

A multivariate analysis, implemented as a boosted decision tree (BDT), is employed to enhance the signal relative to the continuum background. This BDT is based on the $15$ variables described in Table 1.

The discriminating variables can be classified into three groups: (a) $B$ meson variables, related to the reconstruction of the decay vertex and to the collinearity between $\overrightarrow{p}^{B}$ and the flight vector between the production and decay vertices $\overrightarrow{\Delta x}$; (b) variables describing the muons that form the $B$ meson candidate; and (c) variables related to the rest of the event. The selection of the variables aims to maximise the discrimination power of the classifier without introducing significant dependence on the invariant mass of the muon pair.

The same discriminating variables were used in the previous analysis based on the full Run 1 dataset [15]. The removal of individual variables was explored to simplify the BDT input, however, this results inevitably in a significant reduction of the BDT separation power. To minimise the dependence of the classifier on the effects of pile-up, the additional tracks considered to compute the variables $I_{0.7}$, DOCAxtrk and $N^{\mathrm{close}}_{\mathrm{xtrk}}$ are required to be compatible with the primary vertex matched to the dimuon candidate.

The correlations among the discriminating variables were studied in the MC samples for signal and continuum background discussed in Section 2, and in data from the sidebands of the $\mu^{+}\mu^{-}$ invariant mass distribution. There are significant linear correlations among the variables $\chi^{2}_{{\mathrm{PV,DV}}\,xy}$, $L_{xy}$, $|d_{0}|^{\mathrm{max}}$-sig., $|d_{0}|^{\mathrm{min}}$-sig. and $\chi^{2}_{\mu,\mathrm{xPV}}$. The variables IP${}_{B}^{3\mathrm{D}}$, DOCAμμ and $I_{0.7}$ have negligible correlation with any of the others used in the classifier.

The simulated signal sample and the data from the dimuon invariant mass sideband regions are used for training and testing the classifier. As discussed in Section 2, simulated signal samples are corrected for muon reconstruction efficiency differencies between simulation and data, and reweighted according to the distributions of $p_{\text{T}}$ and $|\eta|$ of the dimuon and of the pile-up observed in data. The BDT training is done using the TMVA toolkit [31].

Sideband data are used for the BDT training and optimisation. The sample is subdivided into three randomly selected separate and equally populated subsamples used in turn to train and validate the selection efficiency of three independent BDTs. The resulting BDTs are found to produce results that are statistically compatible, and are combined in one single classifier in such a way that each BDT is applied only to the part of the data sample not involved in the BDT training.

Figure 2 shows the distribution of the BDT output variable for simulated signal and backgrounds, separately for continuum background and partially reconstructed events.

Also shown is the BDT distribution for dimuon candidates from the sidebands of the invariant mass distribution in data. The BDT output was found to not have any significant correlation with the dimuon invariant mass. The final selection requires a BDT output value larger than $0.1439$, corresponding to signal and continuum background efficiencies of 72% and 0.3% respectively. The analysis uses all candidates after this selection; however, accepted events with BDT values close to the selection threshold are effectively only constraining the background models.444The $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ signal fit was found to be insensitive to the signal for candidates with BDT output value smaller than 0.2455 (corresponding to 54% and 0.03% efficiencies for signal and background respectively). For this reason, signal and reference channel yields and efficiencies are measured relative to the signal reference selection discussed in Section 9, while the events in the final selection with lower BDT values are used to improve the background modelling.

7 Data–simulation comparisons

Despite the lack of correlation between the BDT variable and the candidates invariant mass, the change in relative contribution of different background components as a function of the BDT output can produce a variation of the background shape in separate BDT bins. The shape parameters and normalisations of the backgrounds are for this reason determined purely from data: the $b\bar{b}\to\mu^{+}\mu^{-}X$ MC simulation is not used for BDT training, computation of efficiencies or normalisation purposes. The only role of the $b\bar{b}\to\mu^{+}\mu^{-}X$ MC is the validation of the functional forms employed to parameterise the backgrounds on data, as discussed in Section 10.1. The dependence of the background parameterisation on the BDT output variable was investigated in data and shown to be consistent with simulation, even for lower BDT values where prompt backgrounds contribute significantly. This observation is supported by the fair agreement of the BDT input variables distribution between the simulated $b\bar{b}\to\mu^{+}\mu^{-}X$ background and what observed on data sidebands. Figure 3 compares the distributions of two of the most discriminating variables in the continuum background MC sample with data in the dimuon mass sidebands. Agreement with the sideband data is fair.

The distributions of the discriminating variables are also used to compare simulation and data in the $B^{+}\to J/\psi\,K^{+}$ and $B^{0}_{s}\to J/\psi\,\phi$ samples. To perform these comparisons, for each variable the contribution of the background is subtracted from the $B^{+}\to J/\psi\,K^{+}$ ($B^{0}_{s}\to J/\psi\,\phi$) signal. For this purpose, a maximum-likelihood fit is performed to the invariant mass distribution, separately in bins of rapidity and transverse momentum. The fit model used is simpler than the one employed for the extraction of the $B^{+}$ signal for normalisation as described in Section 8, but is sufficient for the purpose discussed here.

Figure 4 shows examples of the distributions of the discriminating variables obtained from data and simulation for the reference samples. Observed differences are used to estimate systematic uncertainties on the efficiency ratio $R_{\varepsilon}=\varepsilon(B^{+}\to J/\psi\,K^{+})/\varepsilon(B^{0}_{(s)}\to\mu^{+}\mu^{-})$ with the procedure described in Section 9. The discrepancy visible for the isolation variable $I_{0.7}$ in the $B^{+}\to J/\psi\,K^{+}$ channel is the most significant among all variables and both reference channels.

8 $B^{+}\rightarrow J/\psi K^{+}$ yield extraction

The reference channel yield is extracted with an unbinned extended maximum-likelihood fit to the $J/\psi K^{+}$ invariant mass distribution. The functional forms used to model both the signal and the backgrounds are obtained from studies of MC samples. All the yields are extracted from the fit to data, while the shape parameters are determined from a simultaneous fit to data and MC samples. Free parameters are introduced for the mass scale and mass resolution to accommodate data–MC differences. The best-fit values indicate a negligibly poorer resolution and a mass shift at the level of $2$ MeV.

The fit includes four components: $B^{+}\to J/\psi\,K^{+}$ decays, Cabibbo-suppressed $B^{+}\to J/\psi\,\pi^{+}$ decays in the right tail of the main peak, partially reconstructed $B$ decays (PRD) where one or more of the final-state particles are missing, and the non-resonant background composed mostly of $b\bar{b}\to J/\psi X$ decays. All components other than the last one have shapes constrained by MC simulation as described below, with the data fit including an additional Gaussian convolution to account for possible data–MC discrepancies in mass scale and resolution. The shape of the $B^{+}\to J/\psi\,K^{+}$ mass distribution is parameterised using a Johnson $S_{U}$ function [32, 33]. The final $B^{+}\to J/\psi\,K^{+}$ yield includes the contribution from radiative effects (i.e. where photons are emitted from the $B$ decay products). The $B^{+}\to J/\psi\,\pi^{+}$ decays are modelled by the sum of a Johnson $S_{U}$ function and a Gaussian function, where all parameters except the normalisation are determined from the simulation. The decay modes contributing to the PRD are classified in simulation on the basis of their mass dependence. Each of the three resulting categories contributes to the overall PRD shape with combinations of Fermi–Dirac and exponential functions, contributing differently in the low-mass region. Their shape parameters are determined from simulation. Finally, the non-resonant background is modelled with an exponential function with the shape parameter extracted from the fit. The normalisation of each component is unconstrained in the fit, which is therefore mostly independent of external inputs for the branching fractions. The residual dependence of the PRD model shapes on the relative branching fractions of the contributing decays is considered as a source of systematic uncertainty. The resulting fit, shown in Figure 5, yields $334\,351$ $B^{+}\to J/\psi\,K^{+}$ decays with a statistical uncertainty of 0.3%. The ratio of yields of $B^{+}\to J/\psi\,\pi^{+}$ and $B^{+}\to J/\psi\,K^{+}$ is $\left(3.71\pm 0.09\right)$% (where the uncertainty reported is statistical only), in agreement with the expectation from the world average [29] of $\left(3.84\pm 0.16\right)$%.

Some systematic uncertainties are included by design in the fit. For example, the effect of the limited MC sample size is included by performing a simultaneous fit to data and MC samples. Scaling factors determined in the fit to data account for the differences in mass scale and resolution between data and simulation. Additional systematic uncertainties are evaluated by varying the default fit model described above. They take into account the kinematic differences between data and the MC samples used in the fit, differences in efficiency between $B^{+}$ and $B^{-}$ decays and uncertainties in the relative fractions and shapes of PRD and in the shape of the various fit components. The stability of this large sample fit is verified by repeating the fit with different initial parameter values. In each case, the change relative to the default fit is recorded, symmetrised and used as an estimate of the systematic uncertainty. The main contributions to the systematic uncertainty come from the functional models of the background components, the composition of the PRD and the signal charge asymmetry. The total systematic uncertainty in the $B^{+}$ yield amounts to $4.8\%$.

9 Evaluation of the $B^{+}\to J/\psi\,K^{+}$ to $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ efficiency ratio

The ratio of efficiencies $R_{\varepsilon}=\varepsilon(B^{+}\to J/\psi\,K^{+})/\varepsilon(B^{0}_{(s)}\to\mu^{+}\mu^{-})$ enters the $\mathcal{D}_{\mathrm{ref}}$ term defined in Section 1: $\mathcal{D}_{\mathrm{ref}}=N_{J/\psi K^{+}}/R_{\varepsilon}$. Both channels are measured in the fiducial acceptance for the $B$ meson, defined as $p_{\text{T}}^{B}>8.0$ GeV and $\left|\eta_{B}\right|<2.5$. Correspondingly, $\varepsilon(B^{+}\to J/\psi\,K^{+})$ and $\varepsilon(B^{0}_{(s)}\to\mu^{+}\mu^{-})$ are measured within the $B$ meson fiducial acceptance and include additional final state particles acceptance as well as trigger, reconstruction and selection efficiencies. The final state particles acceptance is defined by the selection placed on the particles in the final state: $|\eta_{\mu}|<2.5$ and $p_{\text{T}}^{\mu}>6.0$ $(4.0)$ GeV for the leading (trailing) muon $p_{\text{T}}$, $p_{\text{T}}^{K}>1.0$ GeV and $\left|\eta_{K}\right|<2.5$ for kaons. The signal reference BDT selection, defined as BDT $>$ 0.2455, has an efficiency of about $54\%$ ($51\%$) in the signal (reference) channel. The overall efficiency ratio $R_{\varepsilon}$ is 0.1176$\,\pm\,$0.0009 (stat.)$\,\pm\,$0.0047 (syst.), with uncertainties determined as described below.

The ratio $R_{\varepsilon}$ is computed using the mean lifetime of $B_{s}^{0}$ [34, 29] in the MC generator. The same efficiency ratios apply to the $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$  decays, within the MC statistical uncertainty of $0.8\%$. The statistical uncertainties in the efficiency ratios come from the finite number of events available for the simulated samples. The systematic uncertainty affecting $R_{\varepsilon}$ comes from five sources.

The first contribution is due to the uncertainties in the data-driven weights introduced in Section 2, and amounts to 0.8%. This term is assessed by creating alternative datasets using correction factors that are randomly sampled in accord with their nominal values and uncertainties. The RMS value of the distribution of $R_{\varepsilon}$ obtained from these datasets is taken as the systematic uncertainty.

A second contribution of 1.0% is related to the muon trigger and reconstruction efficiencies. The effect of the uncertainties in the data-driven efficiencies is evaluated using random sampling, as above.

A 3.2% systematic uncertainty contribution arises from the differences between data and simulation observed in the modelling of the discriminating variables used in the BDT classifier (Table 1). For each of the $15$ variables, the MC samples for $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ and $B^{+}\to J/\psi\,K^{+}$ are reweighted with the ratio of the $B^{+}\to J/\psi\,K^{+}$ event distributions in sideband-subtracted data and the MC simulation. The isolation variable $I_{0.7}$ is computed using charged-particle tracks only, and differences between $B^{+}$ and $B_{s}^{0}$ are expected and were observed in previous studies [26]. Hence for this variable the reweighting procedure for the $B^{0}_{s}\to\mu^{+}\mu^{-}$ MC sample is based on $B^{0}_{s}\to J/\psi\,\phi$ data. For all discriminating variables except $I_{0.7}$, the value of the efficiency ratio is modified by less than 2% by the reweighting procedure and each variation is taken as an independent contribution to the systematic uncertainty in the efficiency ratio. For $I_{0.7}$ the reweighting procedure changes the efficiency ratio by about $6\%$. Because of the significant mis-modelling, the MC samples obtained after reweighting on the distribution of $I_{0.7}$ are taken as a reference, thus correcting the central value of the efficiency ratio. The 1% uncertainty in the $I_{0.7}$ correction is added to the sum in quadrature of the uncertainties assigned to the other discriminating variables. The total uncertainty in the modelling of the discriminating variables is the dominant contribution to the systematic uncertainty in $R_{\varepsilon}$.

A fourth source of systematic uncertainty arises from differences between the $B^{0}_{s}\to\mu^{+}\mu^{-}$ and the $B^{+}\to J/\psi\,K^{+}$ channel related to the reconstruction efficiency of the kaon track and of the $B^{+}$ decay vertex. These uncertainties are mainly due to inaccuracy in the modelling of passive material in the ID. The corresponding systematic uncertainty is estimated by varying the detector model in simulations, which results in changes between 0.4% and 1.5% depending on the $\eta$ range considered. The largest value is used in the full eta range.

Finally, the uncertainty associated with reweighting the simulated events as a function of the pile-up multiplicity distribution contributes 0.6%.

Table 2 summarises these systematic uncertainties.

The efficiency ratio enters in Eq. (1) with the $\mathcal{D}_{\mathrm{ref}}$ term defined in Section 1, multiplied by the number of observed $B^{\pm}$ candidates. The total uncertainty in $\mathcal{D}_{\mathrm{ref}}$ is $\pm 6.3$%.

A correction to the efficiency ratio for $B^{0}_{s}\to\mu^{+}\mu^{-}$ is needed because of the width difference $\Delta\Gamma_{s}$ between the $B_{s}^{0}$ eigenstates. According to the SM, the decay $B^{0}_{s}\to\mu^{+}\mu^{-}$ proceeds mainly through the heavy state $B_{s,\mathrm{H}}$ [1, 16], which has width $\Gamma_{s,\mathrm{H}}=\Gamma_{s}-\Delta\Gamma_{s}/2$, which is $6.6$% smaller than the average $\Gamma_{s}$ [29]. The variation in the value of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ mean lifetime was tested with simulation and found to change the $B_{s}^{0}$ efficiency, and consequently the $B_{s}^{0}$ to $B^{+}$ efficiency ratio, by $+3.3$%. This correction is applied to the central value of $\mathcal{D}_{\mathrm{ref}}$ used in Section 11 for the determination of $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$).555The decay time distribution of $B^{0}_{s}\to\mu^{+}\mu^{-}$ is predicted to be different from the one of $B_{s,\mathrm{H}}$ in scenarios of new physics, with the effect related to the observable $\mathcal{A}^{\,\mu\mu}_{\,\Delta\Gamma}$ [16, 17]. The maximum possible deviation from the SM prediction of $\mathcal{A}^{\,\mu\mu}_{\,\Delta\Gamma}=+1$ is for $\mathcal{A}^{\,\mu\mu}_{\,\Delta\Gamma}=-1$, for which the decay time distribution of $B^{0}_{s}\to\mu^{+}\mu^{-}$ corresponds to the distribution of the $B_{s,\mathrm{L}}$ eigenstate. In the comparison with new-physics predictions, the value of $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) obtained from this analysis should be corrected by $+3.6\%$ or $+7.8\%$ respectively for $\mathcal{A}^{\,\mu\mu}_{\,\Delta\Gamma}=0$ and $-1$.

Due to the small value of $\Delta\Gamma_{d}$, no correction needs to be applied to the $B^{0}\to\mu^{+}\mu^{-}$ decay.

10 Extraction of the signal yield

Dimuon candidates passing the preliminary selection and the selections against hadron misidentification and continuum background are classified according to four intervals (with boundaries at 0.1439, 0.2455, 0.3312, 0.4163 and 1) in the BDT output. Repeating the Run 1 analysis approach, each interval is chosen to give an equal efficiency of $18\%$ for signal MC events, and they are ordered according to increasing signal-to-background ratio.

An unbinned extended maximum-likelihood fit is performed on the dimuon invariant mass distribution simultaneously across the four BDT intervals, each including models for the respective signal and background contributions. The first two bins (covering the lowest BDT values considered) contribute mostly to background modelling: it has been verified with MC pseudo-experiments that they have negligible impact on the signal extraction.

The result of the fit is parameterised with the total yield of $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$ events in the three highest intervals of BDT output. Section 10.1 describes the signal and background fit models. The parameters describing the background are allowed to vary freely and are determined by the fit. The normalisations of the individual fit components, including the signals, are completely unconstrained and allowed to take negative values. The ratios of the signal yields in different BDT bins are constrained to equal the ratios of the signal efficiencies in those same bins. The systematic uncertainties due to variations in the relative signal and background efficiencies between BDT intervals, to the signal parameterisation and to the background model are discussed in Sections 10.2 and 10.3. Each is modelled in the likelihood as a multiplicative Gaussian distribution whose width is equal to the corresponding systematic uncertainty.

10.1 Signal and background model

The signal and background models are derived from simulations and from data collected in the mass sidebands of the search region.

The invariant mass distribution of the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ signal is described with two double-Gaussian distributions, centred respectively at the $B^{0}$ or $B_{s}^{0}$ mass. The shape parameters are extracted from simulation, where they are found to be uncorrelated with the BDT output. Systematic uncertainties in the mass scale and resolutions are considered separately. Figure 6 shows the invariant mass distributions for $B^{0}$ and $B_{s}^{0}$, obtained from MC events and normalised to the SM expectations.

The background in the signal fit is composed of the types of events described in Section 4: (a) the continuum background; (b) the background from partially reconstructed $b\to\mu^{+}\mu^{-}X$ events, which is present mainly in the low mass sideband; (c) the peaking background.

The non-peaking contributions have a common mass shape model, with parameters constrained across the BDT bins in the fit as described below, and independent yields across BDT bins and components. Systematic uncertainties arising from model assumptions will be discussed in Section 10.3, including effects due to the presence of $B_{c}^{+}\to J/\psi\mu^{+}\nu$ and semileptonic $B^{0}_{(s)}/\Lambda_{b}^{0}\to h\mu\nu$ decays.

Both in simulation and sideband data, the continuum background has a small linear dependence on the dimuon invariant mass. In the simulation, the slope parameter has a roughly linear dependence versus BDT interval; the mass sidebands in data confirm this trend, albeit with large statistical uncertainty. This dependence is included in the fit model. The small systematic uncertainties due to deviations from this assumption are discussed below in Section 10.3.

The $b\to\mu^{+}\mu^{-}X$ background has a dimuon invariant mass distribution that falls monotonically with increasing dimuon mass. The mass dependence is derived from data in the low mass sideband, and described with an exponential function with the same shape in each BDT interval. The value of the shape parameter is extracted from the fit to data.

The invariant mass distribution of the peaking background is very similar to the $B^{0}$ signal, as shown in Figure 1(b). The description of this component is obtained from MC simulation, which indicates that the shape and normalisation are the same for all BDT bins. In the fit, this contribution is included with fixed mass shape and with a normalisation of $2.9\pm 2.0$ events, as discussed in Section 5. This contribution is equally distributed among the three highest intervals of the BDT output.

The fitting procedure is tested with MC pseudo-experiments, as discussed in Section 10.3.

10.2 Relative signal efficiency between BDT bins

Section 9 explains how systematic uncertainties affect the overall selection efficiency for signal candidates. The separation of the candidates according to BDT bins introduces an additional dependence on the relative efficiencies in each BDT bin, and systematic uncertainties in these relative efficiencies must be accounted for. Two different procedures are explored. First, the distribution of the BDT output is compared between MC simulation and background-subtracted data for the reference and control channels. The differences observed in the ratio of data to simulation are described with a linear dependence on the BDT output. The linear dependencies observed for $B^{+}\to J/\psi\,K^{+}$ and $B^{0}_{s}\to J/\psi\,\phi$ are in turn used to reweight the BDT-output distribution in the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ MC sample. The maximum corresponding absolute variations in the efficiencies are equal to $+1.7$% and $-2.3$% respectively in the second and fourth BDT intervals, with the third interval basically unaffected. A second assessment of the systematic uncertainties in the relative efficiency of the BDT intervals is obtained with a procedure similar to the one used for the event selection (Section 9). For each discriminating variable, the MC sample is reweighted according to the difference between simulation and data observed in the reference channels. The variation in the efficiency of each BDT interval is taken as the contribution to the systematic uncertainty due to mis-modelling of that variable. The sum in quadrature of the variations due to all discriminating variables is found to be similar in the $B^{+}\to J/\psi\,K^{+}$ and $B^{0}_{s}\to J/\psi\,\phi$ channels. Absolute variations of $\pm 1.0$%, $\pm 2.4$% and $\pm 4.4$% are found in the second, third and fourth BDT intervals respectively. The first of these procedures is used as a baseline for inclusion of Gaussian terms in the signal extraction likelihood to account for the uncertainty in the relative signal efficiency in the three most sensitive BDT bins. Care is taken in constraining the sum of the efficiencies of the three intervals sensitive to the signal, since that absolute efficiency and the corresponding uncertainty are parameterised with the $R_{\varepsilon}$ term.

Figure 7 shows the distribution of the BDT output from data and simulation for the reference channels, after reweighting the MC sample. The MC distribution for $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ events is also shown, illustrating how the linear deviation obtained from the reference channels affects the simulated signal BDT output. When studying these effects, the linear fits to the ratios in Figures 7(a) and 7(b) are performed in the range corresponding to the three BDT bins with the highest signal-to-background ratio, since the remaining bin is insensitive to the signal contribution.

10.3 Systematic uncertainties in the fit

Studies based on MC pseudo-experiments are used to assess the sensitivity of the fit to the input assumptions. Variations in the description of signal and background components are used in the generation of these samples. The corresponding changes in the average numbers, $N_{s}$ and $N_{d}$, of $B_{s}^{0}$ and $B^{0}$ events determined by the fit, run in the nominal configuration, are taken as systematic uncertainties. The size of the variations used in the generation of the MC pseudo-experiments is determined in some cases by known characteristics of the ATLAS detector (reconstructed momentum scale and momentum resolution), in others using MC evaluation (background due to semileptonic $B^{0}_{(s)}/\Lambda_{b}^{0}\to h\mu\nu$ decays and to $B_{c}^{+}\to J/\psi\mu^{+}\nu$), and in others from uncertainties determined from data in the sidebands or from simulation (shapes of the background components and their variation across the BDT intervals).

The MC pseudo-experiments were generated with the normalisation of the continuum and $b\to\mu\mu X$ components obtained from the fit to the data in the sidebands of the invariant mass distribution, and the peaking background from the expectation discussed in Section 5. The signal was generated with different configurations, roughly covering the range between zero and twice the expected SM yield.

For all variations of the assumptions and all configurations of the signal amplitudes the distributions of the differences between fit results and generated values are used to evaluate systematic uncertainties. In addition, distributions obtained from MC pseudo-experiments generated and fitted according to the nominal fit model are used to study systematic biases deriving from the fit procedure. For both signal yields, the bias is smaller than 15% of the fit error, for true values of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ branching fraction above $5\times 10^{-10}$.

The shifts in $N_{s}$ or $N_{d}$ are combined by considering separately the sums in quadrature of the positive and negative shifts and taking the larger as the symmetric systematic uncertainty. The total systematic uncertainty is found to increase with the assumed size of the signal, with a dependence $\sigma_{\mathrm{syst}}^{N_{s}}=3+0.05N_{s}$ and $\sigma_{\mathrm{syst}}^{N_{d}}=2.9+0.05N_{s}+0.05N_{d}$. Most of the shifts observed have opposite sign for $N_{s}$ and $N_{d}$, resulting in a combined correlation coefficient in the systematic uncertainties of $\rho_{\mathrm{syst}}=-0.83$.

The systematic uncertainties discussed in this Section are included in the fit to the $\mu^{+}\mu^{-}$ candidates in data. The fit for the yield of $B_{s}^{0}$ and $B^{0}$ events is modified by including in the likelihood two smearing parameters for $N_{s}$ and $N_{d}$ that are constrained by a two-dimensional Gaussian distribution parameterised by the values of $\sigma_{\mathrm{syst}}^{N_{s}}$, $\sigma_{\mathrm{syst}}^{N_{d}}$ and $\rho_{\mathrm{syst}}$.

10.4 Results of the signal yield extraction

The numbers of background events contained in the signal region ($5166$–$5526$ MeV) are computed from the interpolation of the data observed in the sidebands. This procedure yields $2685\pm 37$, $330\pm 14$, $51\pm 6$ and $7.9\pm 2.6$ events respectively in the four intervals of BDT output. For comparison, the total expected numbers of signal events according to the SM prediction are $91$ and $10$ for $N_{s}$ and $N_{d}$ respectively, equally distributed among the three intervals with the highest signal-to-background ratio.

In those three BDT intervals, in the unblinded signal region, a total of $1951$ events in the full mass range of $4766$–$5966$ MeV are used in the likelihood fit to signal and background. Without applying any bounds on the values of the fitted parameters, the values determined by the fit are $N_{s}=80\pm 22$ and $N_{d}=-12\pm 20$, where the uncertainties correspond to likelihood variations satisfying $-2\,\Delta\ln(\mathcal{L})=1$. The likelihood includes the systematic uncertainties discussed above, but statistical uncertainties largely dominate. The result is consistent with the expectation from simulation. The uncertainties in the result of the fit are discussed in Section 11, where the measured values of the branching fractions are presented.

Figure 8 shows the dimuon invariant mass distributions in the four BDT intervals, together with the projections of the likelihood. A modified Kolmogorov–Smirnov (KS) test [35] is used to estimate the fit quality: the $p$-value is estimated by comparing the maximum of the KS distance across the four histograms of Figure 8 with the distribution of the same quantity from pseudo-experiments generated with the shape resulting from the fit to data. This procedure yields a compatibility probability of 84%.

11 Branching fraction extraction

The branching fractions for the decays $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $B^{0}\to\mu^{+}\mu^{-}$  are extracted from data using a maximum-likelihood fit. The likelihood is obtained from the one used for $N_{s}$ and $N_{d}$ by replacing the fit parameters with the corresponding branching fractions divided by normalisation terms in Eq. (1), and including Gaussian multiplicative factors for the normalisation uncertainties. All results are obtained profiling the fit likelihood with respect to all parameters involved other than the branching fraction(s) of interest.

The normalisation terms include external inputs for the $B^{+}$ branching fraction and the relative hadronisation probability. The branching fraction is obtained from world averages [29] as the product of ${\cal B}(B^{+}\to J/\psi\,K^{+})=(1.010\pm 0.029)\times 10^{-3}$ and ${\cal B}(J/\psi\to\mu^{+}\mu^{-})=(5.961\pm 0.033)\%$. The relative hadronisation probability is equal to one for $B^{0}$, while for $B_{s}^{0}$ it is taken from the latest available HFLAV average [34] $f_{s}/f_{d}=0.256\pm 0.013$, which assumes $f_{u}/f_{d}=1$.

The efficiency-weighted and luminosity-weighted number of events for the reference channel enters in Eq. (1) with the denominator $\mathcal{D}_{\mathrm{ref}}$. The values $\mathcal{D}_{\mathrm{ref}}=(5.69\pm 0.36)\times 10^{6}$ for $B_{s}^{0}$ and $(5.84\pm 0.37)\times 10^{6}$ for $B^{0}$ are obtained from the results of Sections 8 and 9, and include the correction to the $B^{0}_{s}\to\mu^{+}\mu^{-}$ efficiency due to the lifetime difference between $B_{s,\mathrm{H}}$ and $B_{s}^{0}$. The combination of $B^{+}$ branching fraction, hadronisation probabilities and $\mathcal{D}_{\mathrm{ref}}$, i.e. the single-event sensitivity, is equal to $(4.02\pm 0.35)\times 10^{-11}$ for $B^{0}_{s}\to\mu^{+}\mu^{-}$ and $(1.059\pm 0.074)\times 10^{-11}$ for $B^{0}\to\mu^{+}\mu^{-}$.

The values of the branching fractions that maximise the profile likelihood are ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(3.2\pm 0.9\right)\times 10^{-9}$ and ${\cal B}(B^{0}\to\mu^{+}\mu^{-})=\left(-1.3\pm 2.1\right)\times 10^{-10}$, where the uncertainties, both statistical and systematic combined, are estimated using a Gaussian approximation in the neighbourhood of the likelihood function maximum. Figure 9(a) shows likelihood contours for the simultaneous fit to $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) and $\cal B$($B^{0}\to\mu^{+}\mu^{-}$), for values of $-2\Delta\ln\left(\mathcal{L}\right)$ equal to 2.3, 6.2 and 11.8, relative to the maximum of the likelihood.

Table 3 gives a breakdown of the estimated contributions of systematic and statistical uncertainties. The results are dominated by statistical uncertainties, with the most prominent source of systematic uncertainty coming from fit uncertainties, where the largest contributors are the mass scale and $b\to\mu^{+}\mu^{-}X$ background parameterisation.

Given the statistical regime of the analysis, the likelihood contours of Figure 9(a) cannot be immediately translated into contours with the conventional coverage of one, two and three Gaussian standard deviations. Moreover, the contours extend into regions of negative branching fractions, which are unphysical. In order to address these points, a Neyman construction [36] is employed to obtain the 68.3%, 95.5% and 99.7% confidence intervals in the $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) – $\cal B$($B^{0}\to\mu^{+}\mu^{-}$) plane. This construction yields the contours shown in Figure 9(b). The same construction is used to determine the $68.3\%$ confidence interval for $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) with MC pseudo-experiments, and results in

[TABLE]

Statistical and systematic uncertainties (shown separately in the first expression, and combined in the last) are separated by repeating the likelihood fit after setting all systematic uncertainties to zero. The fitted branching fractions (obtained from the unconstrained likelihood maximum) are in all cases inputs to the Neyman construction, which, by design, results in physically allowed values for the resulting branching fractions.

The upper limit at 95% CL on $\cal B$($B^{0}\to\mu^{+}\mu^{-}$), determined with the same Neyman procedure, is

[TABLE]

Using the predicted SM branching fractions from Section 1, the analysis is expected to yield on average a measurement of $\left(3.6^{+1.1}_{-1.0}\right)\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-9}$ for $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) and an upper limit at 95% CL of $7.1\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-10}$ for $\cal B$($B^{0}\to\mu^{+}\mu^{-}$).

12 Combination with the Run 1 result

The likelihood function from the current result is combined with the likelihood function from the Run 1 result [15]. The only common parameters in the combination are the fitted $\cal B$($B^{0}_{(s)}\to\mu^{+}\mu^{-}$) and the combination of external inputs $\mathcal{F}_{\mathrm{ext}}={\cal B}(B^{+}\to J/\psi\,K^{+})\times{\cal B}(J/\psi\to\mu^{+}\mu^{-})\times(f_{u}/f_{s})=(2.35\pm 0.14)\times 10^{-4}$. Except for $\mathcal{F}_{\mathrm{ext}}$, all nuisance parameters are treated as uncorrelated between the two likelihoods, with both likelihoods including their individual parameterisations of systematic uncertainty effects. A negligible change in the results, corresponding to shifts in central values and uncertainties between 1% and 4%, is found when all sources of systematic uncertainty are assumed to be fully correlated.

The maximum of the combined likelihood corresponds to

[TABLE]

Figure 10 shows the likelihood contours for the combined Run 1 and Run 2 result for $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) and $\cal B$($B^{0}\to\mu^{+}\mu^{-}$), for values of $-2\Delta\ln\left(\mathcal{L}\right)$ equal to 2.3, 6.2 and 11.8, relative to the maximum of the likelihood. The contours for the result from 2015–2016 Run 2 data are overlaid for comparison.

When applying the one-dimensional Neyman construction described in Section 11 to this combined likelihood, whose maximum is unconstrained and allowed to access the unphysical (negative) region, the $68.3\%$ confidence interval obtained for $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) is

[TABLE]

The upper limit at 95% CL on $\cal B$($B^{0}\to\mu^{+}\mu^{-}$) is determined with the same Neyman procedure, yielding

[TABLE]

Using the predicted SM branching fractions from Section 1, the analysis is expected to yield on average a measurement of $\left(3.6^{+0.9}_{-0.8}\right)\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-9}$ for $\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) and an upper limit of $5.6\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-10}$ for $\cal B$($B^{0}\to\mu^{+}\mu^{-}$).

The Run 1 and Run 2 results are found to be 1.2 standard deviations apart. Using both runs, the combined significance of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ signal is estimated to be 4.6 standard deviations, and the combined branching fraction measurements differ by 2.4 standard deviations from the SM values in the $\cal B$($B^{0}\to\mu^{+}\mu^{-}$)–$\cal B$($B^{0}_{s}\to\mu^{+}\mu^{-}$) plane. These significancies are assessed purely from the evaluation of likelihood ratios.

13 Conclusions

A study of the rare decays of $B_{s}^{0}$ and $B^{0}$ mesons into oppositely charged muon pairs is presented, based on $36.2$ fb*-1* of $13$ TeV LHC proton–proton collision data collected by the ATLAS experiment in 2015 and 2016.

For the $B_{s}^{0}$ the branching fraction is determined to be ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(3.2^{+1.1}_{-1.0}\right)\times 10^{-9}$, where the uncertainty includes both the statistical and systematic contributions. The result is consistent with the analysis expectation of $\left(3.6^{+1.1}_{-1.0}\right)\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-9}$ in the SM hypothesis. For the $B^{0}$ an upper limit ${\cal B}(B^{0}\to\mu^{+}\mu^{-})<4.3\times 10^{-10}$ is placed at the 95% confidence level, with an expected upper bound of $7.1\leavevmode\nobreak\ \times\leavevmode\nobreak\ 10^{-10}$ in the SM hypothesis. The limit is compatible with the SM prediction.

The result presented in this paper is combined with the ATLAS result from the full Run 1 dataset to obtain ${\cal B}(B^{0}_{s}\to\mu^{+}\mu^{-})=\left(2.8^{+0.8}_{-0.7}\right)\times 10^{-9}$ and ${\cal B}(B^{0}\to\mu^{+}\mu^{-})<2.1\times 10^{-10}$. All the results presented are compatible with the branching fractions predicted by the SM as well as with currently available experimental results.

Acknowledgements

We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.

We acknowledge the support of ANPCyT, Argentina; YerPhI, Armenia; ARC, Australia; BMWFW and FWF, Austria; ANAS, Azerbaijan; SSTC, Belarus; CNPq and FAPESP, Brazil; NSERC, NRC and CFI, Canada; CERN; CONICYT, Chile; CAS, MOST and NSFC, China; COLCIENCIAS, Colombia; MSMT CR, MPO CR and VSC CR, Czech Republic; DNRF and DNSRC, Denmark; IN2P3-CNRS, CEA-DRF/IRFU, France; SRNSFG, Georgia; BMBF, HGF, and MPG, Germany; GSRT, Greece; RGC, Hong Kong SAR, China; ISF and Benoziyo Center, Israel; INFN, Italy; MEXT and JSPS, Japan; CNRST, Morocco; NWO, Netherlands; RCN, Norway; MNiSW and NCN, Poland; FCT, Portugal; MNE/IFA, Romania; MES of Russia and NRC KI, Russian Federation; JINR; MESTD, Serbia; MSSR, Slovakia; ARRS and MIZŠ, Slovenia; DST/NRF, South Africa; MINECO, Spain; SRC and Wallenberg Foundation, Sweden; SERI, SNSF and Cantons of Bern and Geneva, Switzerland; MOST, Taiwan; TAEK, Turkey; STFC, United Kingdom; DOE and NSF, United States of America. In addition, individual groups and members have received support from BCKDF, CANARIE, CRC and Compute Canada, Canada; COST, ERC, ERDF, Horizon 2020, and Marie Skłodowska-Curie Actions, European Union; Investissements d’ Avenir Labex and Idex, ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; CERCA Programme Generalitat de Catalunya, Spain; The Royal Society and Leverhulme Trust, United Kingdom.

The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada), NDGF (Denmark, Norway, Sweden), CC-IN2P3 (France), KIT/GridKA (Germany), INFN-CNAF (Italy), NL-T1 (Netherlands), PIC (Spain), ASGC (Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [37].