Search for scalar resonances decaying into $\mu^{+}\mu^{-}$ in events with and without $b$-tagged jets produced in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

TL;DR

This paper reports a search for narrow scalar resonances decaying into muon pairs in proton-proton collisions at 13 TeV using ATLAS data, setting upper limits on production cross sections without observing significant signals.

Contribution

First search for scalar resonances decaying into muon pairs in events with and without b-tagged jets at 13 TeV with ATLAS, providing new upper limits.

Findings

No significant excess observed above Standard Model background.

Upper limits on cross section times branching ratio range from 1.6 to 44 fb.

Results are consistent with Standard Model expectations.

Abstract

A search for a narrow scalar resonance decaying into an opposite-sign muon pair produced in events with and without $b$-tagged jets is presented in this paper. The search uses 36.1 fb$^{-1}$ of $\sqrt{s}= 13$ TeV proton-proton collision data recorded by the ATLAS experiment at the LHC. No significant excess of events above the expected Standard Model background is observed in the investigated mass range of 0.2 to 1.0 TeV. The observed upper limits at 95$\%$ confidence level on the cross section times branching ratio for $b$-quark associated production and gluon-gluon fusion are between 1.9 and 41 fb and 1.6 and 44 fb respectively, which is consistent with expectations.

Tables

Table 1. Table 1: Normalization factors of the dominant backgrounds as measured in a fit to data under the background plus signal hypothesis (480 GeV b ​ b ​ Φ 𝑏 𝑏 Φ bb\Phi signal). The uncertainty includes both the statistical and systematic sources: the latter ones dominate. The reduction of the observed relative uncertainty of the Z ​ + HF 𝑍 + HF Z\text{+\,HF} normalization in the data fit is due to the large increase in the number of Z ​ + HF 𝑍 + HF Z\text{+\,HF} events when the data in signal regions are included.

Background	Asimov data	Observed data
$Z\text{+\hspace{0.17em}HF}$	$1.00\pm 0.23$	$1.47\pm 0.26$
$Z\text{+\hspace{0.17em}LF}$	$1.00\pm 0.02$	$1.02\pm 0.02$
$t\overline{t}$	$1.00\pm 0.04$	$1.02\pm 0.04$

Full text

\AtlasTitle

Search for scalar resonances decaying into $\mu^{+}\mu^{-}$ in events with and without $b\text{-tagged}$ jets produced in proton–proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector \AtlasAbstract A search for a narrow scalar resonance decaying into an opposite-sign muon pair produced in events with and without $b\text{-tagged}$ jets is presented in this paper. The search uses 36.1 fb*-1* of $\sqrt{s}$ = 13 TeV proton–proton collision data recorded by the ATLAS experiment at the LHC. No significant excess of events above the expected Standard Model background is observed in the investigated mass range of 0.2 to 1.0 TeV. The observed upper limits at 95% confidence level on the cross section times branching ratio for $b\text{-quark}$ associated production and gluon–gluon fusion are between 1.9 and 41 fb and 1.6 and 44 fb respectively, which is consistent with expectations.

\AtlasRefCodeHIGG-2017-10 \PreprintIdNumberCERN-EP-2018-300 \AtlasJournalRefJHEP 07 (2019) 117 \AtlasDOI10.1007/JHEP07(2019)117 \AtlasCoverSupportingNoteSearch for high mass dimuon signatures in b-tagged and b-veto eventshttps://cds.cern.ch/record/2270631 \AtlasCoverAnalysisTeamSinead Farrington, Kostas Nikolopoulos, Christopher McNicol, Elisabetta Pianori, Paul Thompson, Ligang Xia \AtlasCoverEdBoardMemberJason Nielsen(chair), Vadim Kostyukhin, Paolo Mastrandrea \AtlasCoverEgroupEditorsatlas-HIGG-2017-10-editors@cern.ch \AtlasCoverEgroupEdBoardatlas-HIGG-2017-10-editorial-board@cern.ch

1 Introduction

This paper presents a search, in proton–proton ($pp$) collision data at a centre-of-mass energy of 13 TeV, targeting a scalar particle $\Phi$ with a mass in the range 0.2–1.0 TeV decaying into two opposite-sign muons. The natural width of the particle is assumed to be much narrower than the experimental resolution, such that it is the latter that dominates the distribution of the dimuon invariant mass in signal events. To maximize the sensitivity to the presence of $b\text{-quarks}$, the analysis is performed in two data categories. The first search category requires at least one jet tagged as containing $b$-hadrons ($b\text{-tagged}$) in the final state. The second category requires exactly zero $b$-tagged jets in the event. Simultaneous fits to these two categories are used to search for both the gluon–gluon fusion and $b$-associated production mechanisms of the scalar particle. This work is primarily motivated as a signature-based search; the event selection is not designed to target any specific model and instead facilitates the comparison of the data with predictions from various theoretical models.

The discovery of the 125 GeV Higgs boson at the Large Hadron Collider (LHC) [1, 2] was a significant step in understanding the mechanism of electroweak symmetry breaking (EWSB). Measurements of the properties of this particle [3, 4, 5, 6, 7] are so far consistent with expectations for the Standard Model (SM) Higgs boson [8, 9, 10, 11, 12, 13]. However, the observed Higgs boson could be only a part of an extended scalar sector, as predicted by several theoretical models beyond the SM. For example, such an extended scalar sector is predicted by a class of extensions of the SM known as the two-Higgs-doublet models (2HDM) [14, 15].

The 2HDM posit the existence of three neutral bosons with properties and coupling strengths differing from those of the SM Higgs particle. One of the higher-mass neutral Higgs bosons may couple to muons at a higher rate than to $\tau$-leptons, in contradiction to the SM Yukawa ordering. This is the case in the Flavourful Higgs model [16], for example. The search described in this paper sets limits on the cross-section times branching ratio of $\Phi$ decaying into muon pairs which are more stringent than the equivalent limits for $\Phi$ decaying into $\tau$-lepton pairs [17]. This is due to higher identification efficiency and lower background rate for muons than for $\tau$-leptons, and the invariant-mass resolution being better for muon pairs than for $\tau$-lepton pairs. Similarly, the coupling of these higher-mass neutral Higgs bosons to $b\text{-quarks}$ could be enhanced relative to the SM Higgs boson coupling. Hence, the production of $\Phi$ in association with $b\text{-quarks}$ ($bb\Phi$), shown in Figure 1(a), could be more abundant than the production of $\Phi$ by gluon–gluon fusion (ggF), shown in Figure 1(b). The $\Phi\rightarrow b\bar{b}$ decay mode, where $\Phi$ is produced in associations with $b\text{-quarks}$, is investigated by CMS in Ref. [18].

Additional interest for the search for an excess of events with respect to the SM predictions in the dimuon plus $b\text{-tagged}$ jets final state arises from $Z^{\prime}$ and leptoquark models [19], which could result in either resonant or non-resonant contributions to the $pp\to t\overline{t}\mu^{+}\mu^{-}\to WbW\bar{b}\mu^{+}\mu^{-}$ final state. This work complements the search for $Z^{\prime}$ candidates decaying into muon pairs presented in Ref. [20]: by classifying events in terms of the presence or absence of a $b\text{-tagged}$ jet, it may be sensitive to signatures not observable with an inclusive selection.

Previous similar searches focused explicitly on beyond-the-SM Higgs bosons decaying into muon pairs, in the mass range 90–500 GeV [21, 22] using $\sqrt{s}$ = 7 and 8 TeV data. This work extends the exclusion limits for new scalar resonances with masses up to 1.0 TeV.

The analysis described in this paper makes use of a fit to the observed dimuon invariant mass ($m_{\mu\mu}$) spectra. The presence of a signal produced either in association with $b\text{-quarks}$ or via gluon–gluon fusion is tested with simultaneous fits to the two data regions (with and without a $b\text{-tagged}$ jet) separately for $bb\Phi$ and ggF production modes. No assumptions are made about the relative contributions of the $bb\Phi$ and ggF cross-sections.

This paper is structured as follows. Section 2 describes the ATLAS detector. Section 3 discusses the data and the simulated event samples used to model the signal and the background processes. The event reconstruction is discussed in Section 4, while Section 5 describes the background estimate and introduces the signal interpolation procedure used to model the $m_{\mu\mu}$ distribution for resonance mass hypotheses for which no simulated sample was generated. The event yields and the description of the statistical analysis are discussed in Section 6, followed by Section 7, which provides an overview of the systematic uncertainties. Section 8 summarizes the results.

2 ATLAS detector

ATLAS [23], [24] is a multipurpose particle detector with a forward–backward symmetric cylindrical geometry and near 4$\pi$ coverage in solid angle.111ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the $z$-axis along the beam pipe. The $x$-axis points from the IP to the centre of the LHC ring, and the $y$-axis points upward. Cylindrical coordinates $(r,\phi)$ are used in the transverse plane, $\phi$ being the azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. The angular distance is measured in units of $\Delta R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer.

The inner tracking detector (ID) covers the pseudorapidity range $|\eta|<2.5$, and is surrounded by a superconducting solenoid providing a $2\text{\,}\mathrm{T}$ magnetic field. At small radii, a high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track. It is followed by the silicon microstrip tracker, which provides eight measurement points per track. The silicon detectors are complemented by a gas-filled straw-tube transition radiation tracker, which extends the tracking capability radially with typically 35 measurements per track for particles at $|\eta|=2.0$.

Electromagnetic (EM) calorimetry, within the region $|\eta|<3.2$, is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) sampling calorimeters, with an additional thin LAr presampler covering $|\eta|<1.8$ to correct for energy loss in the upstream material. For $|\eta|<2.5$ the EM calorimeter is divided into three layers in depth. A steel/scintillator-tile calorimeter provides hadronic calorimetry for $|\eta|<1.7$. LAr ionization, with copper as absorber, is used for the hadronic calorimeters in the endcap region $1.5<|\eta|<3.2$. The solid-angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules in $3.1<|\eta|<4.9$, optimized for EM and hadronic measurements, respectively.

The muon spectrometer (MS) surrounds the calorimeters and comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field provided by one barrel and two end-cap toroid magnets. The precision chamber system covers the region $|\eta|<2.7$ with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region. The muon trigger system covers the range $|\eta|<2.4$ with resistive-plate chambers in the barrel and thin-gap chambers in the endcap regions.

A two-level trigger and data acquisition system is used to record events for offline analysis [25]. The level-1 trigger is implemented in hardware and uses a subset of the detector information to reduce the event rate to a value of at most $100\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$. It is followed by a software-based high-level trigger which filters events using the full detector information and outputs events for permanent storage at an average rate of $1\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$.

3 Data and simulated event samples

This search is performed with a sample of $pp$ collision data recorded at a centre-of-mass energy $\sqrt{s}$ =13 TeV corresponding to an integrated luminosity of 36.1 fb*-1* collected during 2015 and 2016. Events are considered for further analysis only if they were collected under stable LHC beam conditions and the relevant detector components were fully operational.

The data used in this analysis were collected using a combination of muon triggers. For the 2015 dataset an isolated muon with transverse momentum $p_{\text{T}}$ greater than $20\,\text{Ge\kern-1.00006ptV}$ was required, while for 2016 this requirement was raised to $26\,\text{Ge\kern-1.00006ptV}$. To avoid inefficiencies due to the isolation requirement, these triggers were complemented by a trigger requiring a reconstructed muon with $p_{\text{T}}$ greater than $50\,\text{Ge\kern-1.00006ptV}$, without any isolation requirements.

Simulated signal events were generated similarly to the SM Higgs boson $H$, where $H$ was replaced by a higher-mass scalar boson $\Phi$ with a constant width of 4 MeV. Nine Monte Carlo (MC) samples were generated in the mass range 0.2–1.0 TeV at intervals of 100 GeV. Both the ggF and $bb\Phi$ production modes were considered. Events were generated at next-to-leading order (NLO) with Powheg-Box 2 [26, 27, 28, 29] and MadGraph5_aMC@NLO [30, 31], for ggF and $bb\Phi$ respectively. The CT10 [32] set of parton distribution functions (PDFs) was used in the generation of ggF events while CT10nlo_nf4 [33] PDFs were used to produce the $b\text{-quark}$ associated signal samples in the four-flavour scheme. In the gluon–gluon fusion production mode, Pythia 8.210 [34] with the AZNLO [35] set of tuned parameters was used together with the CTEQ6L1 [36] PDF set for the parton shower, underlying event and hadronization simulation. For the $b\text{-quark}$ associated production, Pythia 8.210 with the A14 [37] set of tuned parameters was used together with the NNPDF2.3LO [38] PDF set for the parton shower, underlying event, and hadronization simulation. The EvtGen v1.2.0 program [39] was used to model the properties of the bottom and charm hadron decays in all signal samples.

Since the generated signal samples are spaced equally in $m_{\Phi}$ while the experimental resolution increases with $m_{\Phi}$, signal distributions for intermediate mass hypotheses were obtained using an interpolation procedure, described in Section 5.

Events containing $Z/\gamma^{*}\text{+\,jets}$ were simulated using Powheg-Box 2 [28, 40] interfaced to the Pythia 8.186 parton shower model and the CT10 PDF set. The AZNLO tune was used, with PDF set CTEQ6L1 for the modelling of non-perturbative effects. The EvtGen v1.2.0 program was used to model the properties of the bottom and charm hadron decays. Photos++ 3.52 [41] was used for photons radiated from electroweak vertices and charged leptons. Generator-level filters are employed to enhance the fraction of simulated events in the phase-space region that is most relevant for the analysis. Event yields were corrected with a mass-dependent rescaling at next-to-next-to-leading order (NNLO) in the QCD coupling constant, computed with VRAP 0.9 [42] and the CT14NNLO PDF set [33]. Mass-dependent electro-weak (EW) corrections were computed at NLO with mcsanc–1.20 [43], along with photon-induced contributions ($\gamma\gamma\to\ell\ell$ via $t$- and $u$-channel processes) computed with the MRST2004QED PDF set [44].

Additional $Z/\gamma^{*}\text{+\,jets}$ events, used to evaluate the systematic uncertainty from the modelling of the fitted $m_{\mu\mu}$ distribution, were simulated with Sherpa 2.2.1  [45] and with MadGraph5_aMC@NLO interfaced with Pythia 8.186. In Sherpa, $Z/\gamma^{*}\text{+\,jets}$ matrix elements were evaluated with up to two extra partons at NLO, and three or four extra partons were included at LO in QCD. The merging of different parton multiplicities was achieved through a matching scheme based on the CKKW-L [46] [47] merging technique using a merging scale of $Q_{\text{cut}}=20$ GeV. The models used for the parton shower and underlying event are the ones provided internally by Sherpa. The Sherpa 2.2.1 generator adopts a full five-flavour scheme, with massless $b\text{-quarks}$ and $c\text{-quarks}$ in the matrix elements. These quarks are given mass in the final state and massive quarks can be produced directly via the parton shower. The MadGraph5_aMC@NLO v2 generator merges the LO (QCD) matrix-element calculations for $Z/\gamma^{*}\text{+\,jets}$ with up to four additional partons (higher jet multiplicities are modelled by the parton shower algorithm). The merging scheme applied to combine different parton multiplicities is the CKKW-L scheme with a merging scale of $Q_{\text{cut}}=30$ GeV. For the LO matrix-element calculation the NNPDF2.3 LO PDFs were used (with $\alpha_{\text{S}}$ = 0.13). The A14 tune, together with the NNPDF2.3 LO PDFs, was used for the parton shower. Similarly to Sherpa 2.2.1, MadGraph5_aMC@NLO adopts a full five-flavour scheme.

For the generation of $t\bar{t}$ and single top-quark production in the $Wt$- and $t$-channel, the Powheg-Box 2 generator with the CT10 PDF set was used for the matrix element calculations. The parton shower, fragmentation, and the underlying event were simulated using Pythia 6.428 with the CTEQ6L1 PDF set and the corresponding Perugia 2012 tune [48]. The top-quark mass was set to $172.5\,\text{Ge\kern-1.00006ptV}$. For the generation of $t\bar{t}$ events, the $h_{\text{damp}}$ parameter of Powheg-Box, which controls the transverse momentum of the first additional emission beyond the Born configuration, was set to the mass of the top quark. The main effect of this parameter is to regulate the high transverse momentum emission against which the $t\bar{t}$ system recoils. The $t\bar{t}$ production sample was normalized to the predicted production cross-section as calculated with the Top++2.0 program to NNLO in perturbative QCD, including soft-gluon resummation to next-to-next-to-leading-log (NNLL) order [49]. The normalization of the single top-quark event samples used an approximate calculation at NLO in QCD for the $t$-channels [50, 51] while NLO+NNLL predictions were used for the $Wt$-channel [51]. The EvtGen v1.2.0 program was used to describe the properties of the bottom and charm hadron decays.

Diboson processes were modelled using the Sherpa 2.1.1  [45] generator and they were calculated for up to one ($ZZ$) or zero ($WW$, $WZ$) additional partons at NLO and up to three additional partons at leading order (LO) using the Comix [52] and OpenLoops [53] matrix-element generators. The matrix element calculation was merged with the Sherpa internal parton shower [54] using the ME+PS@NLO prescription [55]. The CT10 PDF set was used in conjunction with dedicated parton shower tuning developed by the Sherpa authors. The generator cross-sections, calculated at NLO, were used in this case.

The generated signal samples were simulated with the fast simulation [56, 57] framework of ATLAS, which replaces the full simulation of the electromagnetic and hadronic calorimeters by a parameterized model. The background samples were simulated using the full Geant4-based simulation of the detector [58]. Finally, the simulated events were processed through the same reconstruction software as the data. The effects of pile-up from multiple proton–proton interactions in the same and neighbouring bunch crossings were accounted for by overlaying minimum-bias events simulated using Pythia 8 with the A2 tune [59] and interfaced with the MSTW2008LO PDFs. Simulated events were then reweighted to match the pile-up conditions observed in the data.

4 Event reconstruction

Interaction vertices are reconstructed [60] from tracks measured by the inner detector. The primary vertex is defined as the vertex with the largest $\sum p_{\text{T}}^{2}$ of its associated tracks. It must have at least two associated tracks, each with transverse momentum $p_{\text{T}}>400\,\text{Me\kern-1.00006ptV}$.

Muon candidates are reconstructed from an ID track combined with a track or track segment detected in the muon spectrometer [61]. These two tracks are used as inputs to a combined fit (for $p_{\text{T}}$ less than 300 GeV) or to a statistical combination (for $p_{\text{T}}$ greater than 300 GeV) [61]. The combined fit takes into account the energy loss in the calorimeter and multiple-scattering effects. The statistical combination for high transverse momenta is performed to mitigate the effects of relative ID and MS misalignments. High-$p_{\text{T}}$ muon candidates are required to have at least three hits, one in each of the three layers of precision chambers in the MS. For these muons, specific regions of the MS where the alignment is suboptimal are vetoed as a precaution, as well as the transition region between the MS barrel and endcap ($1.01<|\eta|<1.10$). Further quality requirements are applied to the consistency of the ID and MS momentum measurements, and to the measured momentum uncertainty for the MS track.

The reconstructed muon candidates are required to have transverse momentum, $p_{\text{T}}$, greater than $30\,\text{Ge\kern-1.00006ptV}$ and pseudorapidity $|\eta|<2.5$. They are further required to be consistent with the hypothesis that they originate from the primary vertex by applying selections on the transverse ($d_{0}$) and longitudinal ($z_{0}$) impact parameters, defined relative to the primary vertex position: $|d_{0}/\sigma_{d_{0}}|<3$ and $|z_{0}\sin{\theta}|<0.5$ mm where $\sigma_{d_{0}}$ denotes the uncertainty in the transverse impact parameter. A loose isolation requirement is applied, based on the sum of the momenta of inner detector tracks which lie within a variable-sized cone, with $\Delta R_{\max}=0.3$, around the muon track. This isolation requirement is tuned to yield a 99% efficiency over the full range of muon $p_{\text{T}}$.

Jets are reconstructed from noise-suppressed energy clusters in the calorimeter [62] with the anti-$k_{t}$ algorithm [63, 64] with radius parameter $R=0.4$. The energies of the jets are calibrated using a jet energy scale (JES) correction derived from both simulation and in situ studies using data [65]. All jets are required to have $p_{\text{T}}>25\,\text{Ge\kern-1.00006ptV}$ and $|\eta|<4.5$. A multivariate selection that reduces contamination from pile-up [66] is applied to jets with $p_{\text{T}}<60$ GeV and with $|\eta|<2.4$ utilizing calorimeter and tracking information to separate hard-scatter jets from pile-up jets.

Selected jets in the central region can be tagged as containing $b$-hadrons ($b\text{-tagged}$) by using a multivariate discriminant (MV2c10) [67, 68] that combines information from an impact-parameter-based algorithm, from the explicit reconstruction of a secondary vertex and from a multi-vertex fitter that attempts to reconstruct the full $b$- to $c$-hadron decay chain. At the chosen working point, the algorithm provides nominal light-flavour ($u$,$d$,$s$-quark and gluon) and $c\text{-jet}$ misidentification rates of 3% and 32%, respectively, for an average 85% $b\text{-jet}$ tagging efficiency, as estimated from simulated $t\bar{t}$ events for jets with $p_{\text{T}}>20\,\text{Ge\kern-1.00006ptV}$ and $|\eta|<2.5$. The flavour-tagging efficiencies from simulation are corrected separately for $b$-, $c$- and light-flavour jets, based on the respective data-based calibration analyses [69].

Simulated events were corrected to reflect the muon and jet momentum scales and resolutions, $b\text{-jet}$ identification algorithm calibration, as well as the muon trigger, identification, and isolation efficiencies measured in data.

An overlap removal procedure is applied to avoid a single particle being reconstructed as two different objects as follows. Jets not tagged as $b\text{-jets}$, but which are reconstructed within $\Delta R=0.2$ of a muon, are removed if they have fewer than three associated tracks or if the muon energy constitutes more than 50$\%$ of the jet energy. Muons reconstructed within a cone of size $\Delta R=\min\left(0.4,0.04+10\,\text{Ge\kern-1.00006ptV}/p_{\text{T}}\right)$ around the jet axis of any surviving jet are removed. Jets are also discarded if they are within a cone of size $\Delta R=0.2$ around an electron candidate. Electrons are reconstructed from clusters of energy deposits in the electromagnetic calorimeter that match ID tracks, and are identified as described in Ref. [70]. For electrons with transverse energy $E_{\text{T}}>10\,\text{Ge\kern-1.00006ptV}$ and pseudorapidity of $|\eta|$ < 2.47, a likelihood-based selection is used at the “loose” operating point defined in [70].

Following the overlap removal, jet cleaning criteria are applied to identify jets arising from non-collision sources or noise in the calorimeters, and any event containing such a jet is removed [71].

The missing transverse momentum, $E_{\text{T}}^{\text{miss}}$, is defined as the negative vector sum of the transverse momenta of muons, electrons and jets associated with the primary vertex. A soft term [72] is added to include well-reconstructed tracks matched to the primary vertex that are not associated with any of the objects.

At least two reconstructed muons are required in the event, and the two highest-$p_{\text{T}}$ muons of the event are used to form a dimuon candidate. If these muons do not have opposite charges, the event is rejected.

5 Signal and background estimate

Events satisfying the preselection criteria of Section 4 are considered for further analysis. Depending on the value of the reconstructed dimuon invariant mass and the number of $b\text{-tagged}$ jets, they are classified as being in either a control region, used to measure the rate of the dominant backgrounds ($Z/\gamma^{*}\text{+\,jets}$ and $t\bar{t}$) or a signal region, used to search for the $\Phi\to\mu\mu$ signal. All signal and control regions are designed to be orthogonal.

Selected events are retained in the signal region if the dimuon invariant mass, $m_{\mu\mu}$, exceeds $160$ GeV, while they are retained in the control region if $100\,\text{Ge\kern-1.00006ptV}<m_{\mu\mu}<160\,\text{Ge\kern-1.00006ptV}$. The control regions do not include the $Z$-boson peak; this avoids constraints on systematic uncertainties in a region of phase space which may behave differently with respect to the high mass tails where the signal search is focused.

Events in the signal region are further classified as “SRbTag”, if at least one $b\text{-tagged}$ jet is identified in the event, or “SRbVeto”, otherwise. No further optimization of the signal region definition is considered, in order to minimize the assumptions made about the signal kinematics.

Events in the control region are further classified as “CRbVeto” if there is no $b\text{-tagged}$ jet identified in the event, while remaining control region events are classified as “CRbTag”, if the missing transverse momentum is less than $100\,\text{Ge\kern-1.00006ptV}$, or “CRttbar” otherwise.

The dominant background sources for this signature are muon pairs through $Z/\gamma^{*}\text{+\,jets}$, top-quark pair, and diboson production. $W$+jets and QCD multi-jet events contribute less than 0.01$\%$ to the total background in signal regions, and therefore are neglected. All relevant background contributions were modelled using simulated samples as described in Section 3.

Simulated $Z/\gamma^{*}\text{+\,jets}$ events were categorized depending on the generator-level “truth” labels of the jets in the event. Simulated jets were truth-labelled according to which hadrons with $p_{\text{T}}>5\,\text{Ge\kern-1.00006ptV}$ were found within a cone of size $\Delta R=0.3$ around the jet axis. If a $b$-hadron was found, the jet was truth-labelled as a $b\text{-jet}$. If no $b$-hadron was found, but a $c$-hadron was present, then the jet was truth-labelled as a $c\text{-jet}$. Otherwise the jet was truth-labelled as a light-flavour (i.e., $u$,$d$,$s$-quark, or gluon) jet. The $Z/\gamma^{*}\text{+\,jets}$ events were classified as $Z\text{+\,heavy flavour}$ ($Z\text{+\,HF}$) if a $b\text{-jet}$ or $c\text{-jet}$ was found at generator level and $Z\text{+\,light flavour}$ ($Z\text{+\,LF}$) otherwise. These two categories of simulated events were treated as different background components.

The predictions for the expected numbers of $t\bar{t}$, $Z\text{+\,HF}$ and $Z\text{+\,LF}$ events are adjusted to match the number of events in data in the control regions via a global fit (see Section 6), while the other backgrounds are normalized to their expected cross-section, discussed in Section 3.

The expected invariant mass distribution of the narrow resonance signal $m_{\Phi}$ is modelled with a double-sided Crystal Ball (DSCB) [73] function in the mass range 0.2 $\text{Te\kern-1.00006ptV}<m_{\Phi}<$ 1.0 TeV for all nine simulated MC samples. The width of the $m_{\mu\mu}$ distribution is dominated by the experimental resolution, which can be described by a Gaussian distribution. The power-law asymmetric terms of the DSCB have enough degrees of freedom to model the $m_{\mu\mu}$ tails for signals in the 0.2–1.0 TeV mass range. In order to interpolate the signal parameterization to any mass value in this fit range, second-order polynomial parameterizations of all six signal-shape parameters as a function of $m_{\Phi}$ are obtained from a simultaneous fit to all the generated mass points $m_{\Phi}$, separately for events with and without at least one $b\text{-tagged}$ jet. Since all signal samples are scaled to the same arbitrary cross-section, any differences in the number of events are due to differences in the acceptances of the selection. The fitted normalization is parameterized with a second-order polynomial in a similar way to the other DSCB parameters. As a result of this procedure, the $m_{\mu\mu}$ distribution for an arbitrary hypothesis mass can be found by constructing a DSCB from parameters calculated using the fitted coefficients and the hypothesized mass $m_{\Phi}$. The interpolation is performed separately for each systematic variation, and for the gluon–gluon fusion and $b\text{-quark}$ associated production signal samples. Binned templates are then generated from the DSCB. To validate the interpolation procedure, one simulated sample at a time is removed from the simultaneous fit, and the template generated from the DCSB is compared with the $m_{\mu\mu}$ distribution from the corresponding simulated sample and they agree within 2–5$\%$, the size of the bin-by-bin MC statistical uncertainty. Nine MC samples were generated in the mass range 0.2–1.0 TeV, every 100 GeV. The interpolation procedure was used to generate signal templates every 10 GeV between 0.2 and 0.3 TeV, every 20 GeV between 0.3 and 0.6 TeV, and every 50 GeV between 0.6 and 1.0 TeV.

The acceptance of $b\text{-quark}$ associated production in the SRbTag region varies in the range 11–19$\%$ depending on $m_{\Phi}$. The acceptance in the SRbVeto region, which is inversely correlated with the SRbTag one, varies from 20$\%$ to 15$\%$. For this reason, both SRbTag and SRbVeto are included in the global fit (see Section 6). The acceptance of gluon–gluon fusion in the SRbVeto region varies in the range 31–35$\%$. The acceptance in the SRbTag is less than 2$\%$ for all masses.

6 Statistical analysis

The test statistic $\tilde{q}_{\mu}$ as defined in Ref. [74] is used to determine the probability that the background-only model is compatible with the observed data, to extract the local $p$-value, and, if no hint of a signal is found in this procedure, to derive exclusion intervals using the CL${}_{\text{s}}$ method. The binned likelihood function is built as the product of Poisson probability terms associated with the bins in the $m_{\mu\mu}$ distribution in the 0.16–1.5 TeV range. It depends on the parameter of interest, on the normalization factors of the dominant backgrounds and on additional nuisance parameters (representing the estimates of the systematic uncertainties) that are each constrained by a Gaussian prior. The bin-by-bin statistical uncertainties of the simulated backgrounds are also considered. Separate fits are performed to test the different hypotheses: only $bb\Phi$ signal, only ggF signal, and signal composed of different mixtures of the two production modes. Data are fitted in the 0.16–1.5 TeV range to check that the background description is correct and to constrain systematic uncertainties, while the signal presence is tested in a narrower 0.2–1.0 TeV range.

Logarithmic binning in $m_{\mu\mu}$ is chosen to scale with experimental dimuon mass resolution (defined as full-width at half-maximum), which varies from 5% at $m_{\mu\mu}=200\,\text{Ge\kern-1.00006ptV}$ to 14% at $m_{\mu\mu}=1\,\text{Te\kern-1.00006ptV}$, while ensuring that the number of simulated background events in each bin is sufficient to reliably predict the background. The numbers of events in the CRbVeto, CRbTag, and CRttbar control regions are also included in the maximum-likelihood fit. Data in these regions are used to constrain the normalizations of the dominant background processes: $Z\text{+\,LF}$, $Z\text{+\,HF}$, and $t\bar{t}$. The normalization factors used to scale the expected number of events for each process to match the data are free parameters in the fit. Their expected uncertainty is estimated using an Asimov dataset [74], a pseudo-data distribution equal to the background plus signal expectation for a given value of the signal cross-section times branching ratio. The expected and measured normalization factors are shown in Table 1. They are measured in a fit to data under the background plus signal hypothesis, where the signal corresponds to $bb\Phi$ production with mass $m_{\Phi}$ = 480 GeV.

As discussed in Section 7, the experimental uncertainties are treated as fully correlated among different processes and regions, while the theoretical uncertainties are mostly uncorrelated among processes.

The numbers of observed events in the signal regions together with the predicted event yields from signal and background processes are shown in Table  2. The numbers are the results of the maximum-likelihood fit to data under the background plus signal hypothesis, where the signal corresponds to $bb\Phi$ production with $m_{\Phi}$ = 480 GeV. The number of predicted signal events corresponds to the expected upper limit (UL) times the cross-section: this highlights the comparison between signal and background at the edge of the analysis sensitivity.

The observed numbers of events are compatible with those expected from SM processes, within uncertainties. The $m_{\mu\mu}$ distribution in the SRbTag region is shown in Figure 2(a) for a $bb\Phi$-only fit. The $m_{\mu\mu}$ distribution in the SRbVeto region is shown in Figure 2(b), for a ggF-only fit. Each background process is normalized according to its post-fit cross-section. The templates for the $m_{\Phi}$ =200 GeV, $m_{\Phi}$ =480 GeV and $m_{\Phi}$ =1 TeV mass hypotheses are normalized to the expected upper limit.

7 Systematic uncertainties

The sources of systematic uncertainty can be divided into three groups: those of experimental origin, those related to the modelling of the simulated backgrounds, and those associated with signal simulation. The finite size of the simulated background samples is also an important source of uncertainty.

7.1 Experimental uncertainties

The dominant experimental uncertainties originate from residual mismodelling of the muon reconstruction and selection after the simulation-to-data efficiency correction factors have been applied. They include the uncertainty obtained from $Z\to\mu\mu$ data studies and a high-$p_{\text{T}}$ extrapolation uncertainty corresponding to the decrease in the muon reconstruction and selection efficiency with increasing $p_{\text{T}}$ which is predicted by the MC simulation. The degradation of the muon reconstruction efficiency was found to be approximately $3\%$ per TeV as a function of muon $p_{\text{T}}$. Uncertainties in the isolation and trigger efficiencies of muons [61], along with the uncertainty in their energy scale and resolution, are estimated from $Z\to\mu\mu$ data taken at $\sqrt{s}$ =13 TeV. These are found to have only a small impact (they account for $<0.4\%$ of the total uncertainty on the fitted value of the signal cross-sections). Other sources of experimental uncertainties are the residual differences in the modelling of the flavour-tagging algorithm, the jet energy scale and the jet energy resolution after the simulation-to-data correction factors are applied. Flavour-tagging simulation-to-data efficiency correction factors are derived [67] separately for $b\text{-jets}$, $c\text{-jets}$, and light-flavour jets. All three correction factors depend on jet $p_{\text{T}}$ (or $p_{\text{T}}$ and $|\eta|$) and are affected by uncertainties from multiple sources. These are decomposed into uncorrelated components which are then treated independently, resulting in three uncertainties for $b\text{-jets}$ and for $c\text{-jets}$, and five for light-flavour jets. The approximate size of the uncertainty in the tagging efficiency is 2% for $b\text{-jets}$, 10% for $c\text{-jets}$ and 30% for light-flavour jets. Additional uncertainties are considered in the extrapolation of the $b\text{-jet}$ efficiency calibration above $p_{\text{T}}=300\,\text{Ge\kern-1.00006ptV}$ and in the misidentification of hadronically decaying $\tau$-leptons as $b\text{-jets}$. The uncertainties in the jet energy scale and resolution are based on their respective measurements in data [65]. The many sources of uncertainty in the jet energy scale correction are decomposed into 21 uncorrelated components which are treated as being independent of one another. An additional uncertainty that specifically affects the energy calibration of $b$- and $c\text{-jets}$ is considered.

The uncertainties in the energy scale and resolution of the jets and leptons are propagated to the calculation of $E^{\text{miss}}_{\text{T}}$, which also has additional uncertainties from the scale, resolution, and efficiency of the tracks used to define the soft term, along with the modelling of the underlying event.

The pile-up modelling uncertainty is assessed by varying the number of pile-up interactions in simulated events. The variations are designed to cover the uncertainty in the ratio of the predicted and measured cross-section of non-diffractive inelastic events producing a hadronic system of mass $m_{X}>13$ GeV [75].

The uncertainty in the combined 2015+2016 integrated luminosity is 2.1$\%$. It is derived, following a methodology similar to that detailed in Ref. [76], and using the LUCID-2 detector for the baseline luminosity measurements [77], from calibration of the luminosity scale using $x$–$y$ beam-separation scans.

7.2 Theoretical uncertainties

Background processes

As discussed earlier, the rates of major backgrounds are measured using data control regions, thus theoretical uncertainties in the predicted cross-section for $Z/\gamma^{*}\text{+\,jets}$ and $t\bar{t}$ processes are not considered. Instead, the effects of theoretical uncertainties in the modelling of the $m_{\mu\mu}$ distribution and in the ratio of the event yields in signal regions to those in the control regions are evaluated. The dominant theoretical uncertainties are the shape uncertainty in the dimuon mass spectra of the $t\bar{t}$ and $Z\text{+\,HF}$ contributions.

For both the $Z\text{+\,LF}$ and $Z\text{+\,HF}$ processes, the following sources of theoretical and modelling uncertainties are considered: PDF uncertainties (estimated via eigenvector variations and by comparing different PDF sets), limited accuracy of the fixed-order calculation (estimated by QCD scale variations), variations in the choice of strong coupling constant value ($\alpha_{\text{S}}(M_{Z})$), EW corrections, and photon-induced corrections. These variations are treated as fully correlated between $Z\text{+\,LF}$ and $Z\text{+\,HF}$ processes.

The PDF variation uncertainty is obtained using the $90\%$ confidence level (CL) CT14NNLO PDF error set and by following the procedure described in Refs. [78, 79]. Rather than a single nuisance parameter to describe the 28 eigenvectors of this PDF error set, which could lead to an underestimation of its effect, a re-diagonalized set of 7 PDF eigenvectors was used [33]. This represents the minimal set of PDF eigenvectors that maintains the necessary correlations, and the sum in quadrature of these eigenvectors matches the original CT14NNLO error envelope well. They are treated as separate nuisance parameters. The uncertainties due to the variation of PDF scales and $\alpha_{\text{S}}$ are derived using VRAP. The former is obtained by varying the renormalization and factorization scales of the nominal CT14NNLO PDF up and down simultaneously by a factor of two. The value of $\alpha_{\text{S}}$ used (0.118) is varied by $\pm 0.002$. The EW correction uncertainty was assessed by comparing the nominal additive ($1+\delta_{\mathrm{EW}}+\delta_{\mathrm{QCD}}$) treatment with the multiplicative approximation (($1+\delta_{\mathrm{EW}}$)($1+\delta_{\mathrm{QCD}}$)) treatment of the EW correction in the combination of the higher-order EW and QCD effects. The uncertainty in the photon-induced correction is calculated from the uncertainties in the quark masses and the photon PDF. Following the recommendations of the PDF4LHC forum [79], an additional uncertainty due to the choice of nominal PDF set is derived by comparing the central values of CT14NNLO with those from other PDF sets, namely MMHT14 [80] and NNPDF3.0 [81]. The maximum width of the envelope of these comparisons is used as the PDF choice uncertainty, but only if it is larger than the width of the CT14NNLO PDF eigenvector variation envelope.

An additional modelling uncertainty is considered for the $Z\text{+\,HF}$ process. The $m_{\mu\mu}$ spectrum simulated by Powheg was compared with those simulated with Sherpa 2.2.1  [45] and with MadGraph5_aMC@NLO interfaced with Pythia 8.186. A functional form was chosen to describe the envelope of the differences in the $m_{\mu\mu}$ distribution shape between events with at least one $b\text{-quark}$ simulated by Sherpa 2.2.1 and MadGraph5_aMC@NLO v2.

These systematic uncertainties not only imply uncertainties in the shape of the $m_{\mu\mu}$ distribution in the signal region, but they also affect the ratio of the expected number of events in the signal region to the expected number in the control region. The effect on the $Z\text{+\,LF}$ normalization in the bVeto signal region is  $2\%$, while the size of the uncertainty for $Z\text{+\,HF}$ in the bTag signal region is $7\%$.

Moreover, the ratio of the expected number of $Z\text{+\,LF}$ events in the bVeto control region to the events in bTag control and signal regions was estimated with the nominal Powheg+Pythia $Z$+jet sample, and with the alternative MadGraph5_aMC@NLO v2 sample. To cover the observed difference between the two, an additional uncertainty of $27\%$ in the normalization of the $Z\text{+\,LF}$ process in the bTag signal and control region was applied.

For the $t\bar{t}$ process, theoretical uncertainties in the modelling of the $m_{\mu\mu}$ spectrum and in the extrapolation from control region to signal region are also considered. These are estimated by comparing the nominal prediction with alternative ones. To estimate the impact of initial and final state radiation modelling, two alternative Powheg+Pythia samples were generated with the following parameters. In the first one, the renormalization ($\mu_{\textsc{r}}$) and factorization ($\mu_{\textsc{f}}$) scale were varied by a factor of 0.5, the value of $h_{\text{damp}}$ was doubled ($2m_{t}$) and the corresponding Perugia 2012 radiation tune variation was used. In the second sample, the renormalization and factorization scales were varied by a factor of 2, while the $h_{\text{damp}}$ parameter was not changed. The associated variation from the Perugia 2012 tune was used. These choices of parameters have been shown to encompass the cases where $\mu_{\textsc{r}}$ and $\mu_{\textsc{f}}$ are varied independently, and covered the measured uncertainties of the data for unfolded $t\bar{t}$ distributions [82, 83]. Differences in parton shower and hadronization models were investigated using a sample where Powheg-Box was interfaced with Herwig++ 2.7.1 [84] with the UE-EE-5 tune [37] and the corresponding CTEQ6L1 PDFs. The nominal sample was also compared with a sample generated with MadGraph5_aMC@NLO 2.2.1 [30] interfaced with Herwig++. A NLO matrix element and CT10 PDF were used for the $t\bar{t}$ hard-scattering process. The parton shower, hadronization and the underlying events were modelled using the Herwig++ 2.7.1 generator. The UE-EE-5 tune and the corresponding CTEQ6L1 PDF were used. A functional form was chosen to encompass the differences in the $m_{\mu\mu}$ distribution shape between the nominal sample and all these alternative samples. These samples were also used to estimate the 3.5$\%$ extrapolation uncertainty from CRttbar to the SRbTag.

Theoretical uncertainties that arise from higher-order contributions to the cross-sections and PDFs and affect the values of the predicted cross-sections for the diboson and single-top backgrounds are considered.

Signal processes

Uncertainties related to signal modelling include the uncertainties associated with the initial- and final-state radiation, the modelling of underlying events, the choice of the renormalization and factorization scales, and the parton distribution functions. None of them has a sizeable effect on the shape of the $m_{\mu\mu}$ spectra, but some of them do affect the acceptance in the signal regions. Uncertainties for the different mass hypotheses were evaluated, but for simplicity, only the largest of the values is used. In the calculation of the uncertainties, appropriate requirements on the muon and jets kinematics are applied at hadron level for the SRbTag and SRbVeto regions.

The factorization and renormalization scales were varied by a factor of two up and down, including correlated and anti-correlated variations, both in Powheg-Box and in MG5_aMC@NLO. For the bb$\Phi$ process, the largest deviation from the value of the nominal acceptance is considered as the final scale uncertainty (2$\%$ in SRbTag, and 1$\%$ in SRbVeto). A 25$\%$ uncertainty for the acceptance of ggF events in the SRbTag is considered, as well as its anti-correlated effect in the SRbVeto region, following the procedure adopted in Ref. [85].

The PDF uncertainties are estimated by taking the envelope of the changes in the acceptance associated with the PDF4LHC15_nlo_100 (PDF4LC15_nlo_nf4_30)  [79] error set for the ggF (bb$\Phi$) signal process. They correspond to changes in acceptances not larger than 1$\%$.

Systematic variations of the parameters of the A14 (for bb$\Phi$) and AZNLO (for ggF) tunes are used to account for the uncertainties associated with the initial and final state radiation and the modelling of underlying events. For bb$\Phi$, these uncertainties correspond to changes of 3.8$\%$ and 3.2$\%$ of the acceptance in the SRbTag and SRbVeto regions. They are associated with a migration of events from one signal region to the other, hence they are treated as anti-correlated between the two signal regions. For ggF, the uncertainty in the SRbVeto is negligible ($<0.5\%$), while it is 3.8$\%$ in the SRbTag.

As discussed in Section 3, the signal samples were simulated with the fast simulation framework of ATLAS, which replaces the full simulation of the electromagnetic and hadronic calorimeters by a parameterized model. No significant differences in acceptance and dimuon invariant mass were observed between a limited number of fast simulation signal samples and samples generated with full Geant4-based simulation, hence no systematic uncertainty is considered for the use of fast simulation.

The current result is dominated by the statistical uncertainty of the data (which accounts for 66$\%$ of the total uncertainty on the fitted value of the signal cross-section), and the systematic uncertainty due to the finite size of the simulated background samples (which accounts for 25$\%$ of the total uncertainty). Amongst the remaining systematic uncertainties, the modelling of the shape of the $m_{\mu\mu}$ distribution for the $Z$+jets and top–antitop backgrounds dominates (3$\%$ of the total uncertainty) followed by the muon identification efficiency (1.5$\%$). All other considered sources of experimental and theoretical uncertainties (which individually have an impact $<0.4\%$) have a combined contribution of less than 5$\%$ to the total uncertainty on the upper limit.

8 Results

The observed $p$-values as a function of $m_{\Phi}$, obtained applying the statistical analysis described in Section 6, are shown in Figure 3(a) for the $bb\Phi$-only fit and in Figure 3(b) for the ggF-only fit. The lines in each figure correspond to separate fits. The dotted line represents the $p$-value for a fit including only SRbTag and all control regions (CRs), the dashed line represents the $p$-value for a fit including only SRbVeto+CRs, while the solid line corresponds to the combined SRbTag+SRbVeto+CRs fit. The dotted lines in the two figures show similar behaviour and so do the dashed lines, while the combined $p$-values are calculated for separate signal production modes implying very different contributions from the SRbVeto and SRbTag regions, which lead to substantially different limits. The largest excess of events above the expected background is observed for the $b\text{-quark}$ associated production at about $m_{\Phi}=480\,\text{Ge\kern-1.00006ptV}$ and amounts to a local significance of $2.3\sigma$, which becomes $0.6\sigma$ after considering the look-elsewhere effect over the mass range 0.2–1.0 TeV. The look-elsewhere effect is estimated using the method described in Ref. [86].

Since the data are in agreement with the predicted backgrounds the results are given in terms of exclusion limits. These are set using the modified frequentist CL${}_{\textrm{s}}$ method [87]. Upper limits on the cross section times branching ratio for a massive scalar particle are set at the 95% confidence level (CL) as a function of the particle mass. The upper limits at 95% CL on $\sigma_{\Phi}\times\mathcal{B}(\Phi\to\mu\mu)$ (where $\sigma_{\Phi}$ is the cross section and $\mathcal{B}$ is the branching ratio) assume the natural width of the resonance to be negligible compared to the experimental resolution, and they cover the mass range 0.2–1.0 TeV. In Figure 4(a) and 4(b) the 95% CL upper limits are shown for $b\text{-quark}$ associated production and gluon–gluon fusion, respectively. The expected limits are in the ranges 1.3–$25\text{\,}\mathrm{fb}$ and 1.8–$25\text{\,}\mathrm{fb}$, respectively. The observed limits are in the ranges 1.9–$41\text{\,}\mathrm{fb}$ and 1.6–$44\text{\,}\mathrm{fb}$, respectively. Both the expected and observed limits are calculated in the asymptotic approximation [74], which was verified with MC pseudo-experiments for the $m_{\Phi}=1000$ GeV hypothesis. The upper limit on the ggF production mode is dominated by the SRbVeto, while for the $bb\Phi$ production mode both the SRbTag and SRbVeto regions contribute to the result because of the relative acceptance of the $bb\Phi$ production mode in the two regions. Limits shown in Figure 4 are obtained in simultaneous fits of the SRbVeto and SRbTag data regions and correspond to the different signal production modes. Usage of the same data events in both the $bb\Phi$-only and the ggF-only fits explains correlations in the behaviour of the observed limits.

Figures 5(a) and 5(b) show the observed and expected 95% CL upper limits on the production cross section times branching ratio for $\Phi\to\mu\mu$ as a function of the fractional contribution from $b\text{-quark}$ associated production (\sigma_{\text{bb\Phi}}/[\sigma_{\text{{bb\Phi}}}+\sigma_{\text{ggF}}]) and the scalar resonance mass. These fractions are derived assuming that other production mechanisms do not affect the result.

A complete set of tables and figures are available at the Durham HepData repository [88].

9 Conclusion

The ATLAS detector at the LHC has been used to search for a massive scalar resonance decaying into two opposite-sign muons, produced via $b\text{-quark}$ associated production or via gluon–gluon fusion, assuming that the natural width of the resonance is negligible compared to the experimental resolution. The search is conducted with 36.1 fb*-1* of $pp$ collision data at $\sqrt{s}=13$ TeV, recorded during 2015 and 2016. The observed dimuon invariant mass spectrum is consistent with the Standard Model prediction within uncertainties for events both with and without a $b\text{-tagged}$ jet over the 0.2–1.0 TeV range. The observed upper limits at 95% confidence level on the cross-section times branching ratio for $b\text{-quark}$ associated production and gluon–gluon fusion are between 1.9 and 41 fb and 1.6 and 44 fb respectively, which is consistent with expectations.

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. [89].