Search for bottom-squark pair production with the ATLAS detector in final states containing Higgs bosons, $b$-jets and missing transverse momentum

TL;DR

This study searches for bottom-squark pair production in proton-proton collisions, focusing on final states with Higgs bosons, b-jets, and missing transverse momentum, setting new exclusion limits up to 1.5 TeV.

Contribution

It presents the first ATLAS search targeting bottom-squarks decaying via Higgs bosons with no observed excess, extending exclusion limits to 1.5 TeV.

Findings

No significant excess over Standard Model background.

Excluded bottom-squark masses up to 1.5 TeV.

Set 95% confidence level limits on supersymmetric models.

Abstract

The result of a search for the pair production of the lightest supersymmetric partner of the bottom quark ($\tilde{b}_{1}$) using 139 fb$^{-1}$ of proton-proton data collected at $\sqrt{s} = 13$ TeV by the ATLAS detector is reported. In the supersymmetric scenarios considered both of the bottom-squarks decay into a $b$-quark and the second-lightest neutralino, $\tilde{b}_{1} \rightarrow b + \tilde{\chi}^{0}_{2}$. Each $\tilde{\chi}^{0}_{2}$ is assumed to subsequently decay with 100% branching ratio into a Higgs boson ($h$) like the one in the Standard Model and the lightest neutralino: $\tilde{\chi}^{0}_{2} \rightarrow h + \tilde{\chi}^{0}_{1}$. The $\tilde{\chi}^{0}_{1}$ is assumed to be the lightest supersymmetric particle (LSP) and is stable. Two signal mass configurations are targeted: the first has a constant LSP mass of 60 GeV; and the second has a constant mass difference between the $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$ of 130 GeV. The final states considered contain no charged leptons, three or more $b$-jets, and large missing transverse momentum. No significant excess of events over the Standard Model background expectation is observed in any of the signal regions considered. Limits at the 95% confidence level are placed in the supersymmetric models considered, and bottom-squarks with mass up to 1.5 TeV are excluded.

Tables

Table 1. Table 1: Definitions for the SRA, alongside the three varying m eff subscript 𝑚 eff m_{\mathrm{eff}} intervals used. The letter appended to the SRA label corresponds to the low (-L), medium (-M) or high (-H) m eff subscript 𝑚 eff m_{\mathrm{eff}} selection. This selection is sensitive to the bulk regions of both signal scenarios. The jets and b 𝑏 b -jets are ordered by p T subscript 𝑝 T p_{\mathrm{T}} .

Variable	SRA	SRA-L	SRA-M	SRA-H
$N_{\mathrm{leptons}}$ (baseline)	= 0	= 0
$N_{\mathrm{jets}}$	$\ge 6$	$\ge 6$
$N_{b\text{-jets}}$	$\ge 4$	$\ge 4$
$E_{\text{T}}^{\text{miss}}$ [GeV]	$>350$	$>350$
$\min \Delta \varphi ({\mathrm{jet}}_{1-4},{𝐩}_T^{\mathrm{miss}})$ [rad]	$>0.4$	$>0.4$
$\tau$ veto	Yes	Yes
$p_T\left(b_1\right)$ [GeV]	$>200$	$>200$
$\Delta R_{\max }(b,b)$	$>2.5$	$>2.5$
$\Delta R_{\text{max-min}}(b,b)$
$m\left(h_{\mathrm{cand}}\right)$ [GeV]	$>$ 80	$>$ 80
$m_{\mathrm{eff}}$ [TeV]	$>1.0$	$\in [1.0,1.5]$	$\in [1.5,2.0]$	$>2.0$

Table 2. Table 2: Definitions for SRB, targeting the compressed region of the Δ ​ m ​ ( χ ~ 2 0 , χ ~ 1 0 ) = 130 Δ 𝑚 subscript superscript ~ 𝜒 0 2 subscript superscript ~ 𝜒 0 1 130 \Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=130   GeV  scenario. The jets and b 𝑏 b -jets are ordered by p T subscript 𝑝 T p_{\mathrm{T}} .

Variable	SRB
$N_{\mathrm{leptons}}$ (baseline)	= 0
$N_{\mathrm{jets}}$	$\ge 5$
$N_{b\text{-jets}}$	$\ge 4$
$E_{\text{T}}^{\text{miss}}$ [GeV]	$>350$
$\min \Delta \varphi ({\mathrm{jet}}_{1-4},{𝐩}_T^{\mathrm{miss}})$ [rad]	$>$ 0.4
$\tau$ veto	Yes
$m{(h_{\mathrm{cand1}},h_{\mathrm{cand2}})}_{\mathrm{avg}}$ [GeV]	$\in [75,175]$
Leading jet not $b$-tagged	Yes
$p_T\left(j_1\right)$ [GeV]	$>350$
$\left|\Delta \varphi (j_1,E_{\text{T}}^{\text{miss}})\right|$ [rad]	$>$ 2.8
$m_{\mathrm{eff}}$ [TeV]	$>1$

Table 3. Table 3: Definitions for SRC, alongside the four varying 𝒮 𝒮 \mathcal{S}  intervals used. The letter appended to the SRC label corresponds to the lower bound on the 𝒮 𝒮 \mathcal{S}  interval. SRC targets small mass splittings between the b ~ 1 subscript ~ 𝑏 1 \tilde{b}_{1} and χ ~ 2 0 subscript superscript ~ 𝜒 0 2 \tilde{\chi}^{0}_{2} , in the m ​ ( χ ~ 1 0 ) = 60 𝑚 subscript superscript ~ 𝜒 0 1 60 m(\tilde{\chi}^{0}_{1})=60   GeV  signal scenario. The jets and b 𝑏 b -jets are ordered by p T subscript 𝑝 T p_{\mathrm{T}} .

Variable	SRC	SRC22	SRC24	SRC26	SRC28
$N_{\mathrm{leptons}}$ (baseline)	= 0	= 0
$N_{\mathrm{jets}}$	$\ge 4$	$\ge 4$
$N_{b\text{-jets}}$	$\ge 3$	$\ge 3$
$E_{\text{T}}^{\text{miss}}$ [GeV]	$>250$	$>250$
$\min \Delta \varphi ({\mathrm{jet}}_{1-4},{𝐩}_T^{\mathrm{miss}})$ [rad]	$>0.4$	$>0.4$
$𝒮$	$>22$	$\in [22,24]$	$\in [24,26]$	$\in [26,28]$	$>28$

Table 4. Table 4: Summary of all control and validation region definitions used in the analysis.

		Control Regions				Validation Regions
Variable	Units	CRA$1\ell$	CRB$1\ell$	CRC$1\ell$	CRC$2\ell$	VRA$0\ell$	VRB$0\ell$	VRC$0\ell$-T	VRC$0\ell$-Z
$E_{\text{T}}^{\text{miss}}$ Trigger		✓	✓	✓	-	✓	✓	✓	✓
Lepton Trigger		-	-	-	✓	-	-	-	-
$E_{\text{T}}^{\text{miss}}$	[GeV]	$>250$	$>300$	$>250$		$>350$	$>350$	$>250$	$\in [250,600]$
$\min \left[\Delta \varphi ({\mathrm{jet}}_{1-4},E_{\text{T}}^{\text{miss}})\right]$	[rad]	-	-	-	$>0.2$	$>0.4$	$>0.4$	$\in [0.2,0.4]$	$>1.2$
$N_{\mathrm{leptons}}$ (baseline)		= 1	= 1	= 1	= 2	= 0	= 0	= 0	= 0
$N_{\mathrm{leptons}}$ (signal)		= 1	= 1	= 1	= 2(SFOS)	-	-	-	-
$p_T\left({\ell }_1\right)$	[GeV]	$>20$	$>20$	$>20$	$>27$	-	-	-	-
$p_T\left({\ell }_2\right)$	[GeV]	-	-	-	$>20$	-	-	-	-
$m_T$	[GeV]	$>20$	$>20$	$>20$	-	-	-	-	-
$m_{\ell \ell }$	[GeV]	-	-	-	$\in [86,106]$	-	-	-	-
$\tau$ veto		-	✓	-	-	✓	✓	-	-
$N_{\mathrm{jets}}$		$\ge 6$	$\ge 4$	$\ge 4$	$\ge 4$	$\ge 6$	$\ge 4$	$\ge 4$	$\ge 4$
$N_{b\text{-jets}}$		$\ge 4$	$\ge 4$	$\ge 3$	$\ge 3$	$=3$	$=3$	$\ge 3$	$=2$
$p_T\left(b_1\right)$	[GeV]	$>100$	-	-	-	$>100$	-	-	-
$p_T\left(j_1\right)$	[GeV]	-	$>350$	-	-	-	$>350$	-	-
Leading jet not b-tagged		-	✓	-	-	-	✓	-	-
$\left|\Delta \varphi (j_1,E_{\text{T}}^{\text{miss}})\right|$	[rad]	-	$>2.2$	-	-	-	$>2.8$	-	-
${\tilde{E}}_T^{\mathrm{miss}}$	[GeV]	-	-	-	$>250$	-	-	-	-
$𝒮$		-	-	$>17$	-			$\in [15,22]$	$\in [20,22]$
$m_{\mathrm{eff}}$	[TeV]	$>1.0$	$>1.0$	-	-	$>1.0$	$>1.0$	-	-

Table 5. Table 5: Expected background event yields and dominant systematic uncertainties on background estimates in the A-type (inclusive), B-type and C-type (inclusive) regions. Individual uncertainties can be correlated, and do not necessarily add up quadratically to the total background uncertainty. The percentages show the size of the uncertainty relative to the total expected background.

Region	SRA		SRB		SRC
Total background expectation	17.1		3.3		37.9
Total background uncertainty	$2.8$	(16%)	$0.9$	(27%)	$6.2$	(16%)
Systematic, experimental	$1.4$	(8%)	$0.3$	(10%)	$3.0$	(8%)
Systematic, theoretical	$2.3$	(13%)	$0.6$	(18%)	$3.2$	(8%)
Statistical, MC samples	$0.7$	(4%)	$0.4$	(12%)	$2.0$	(5%)

Table 6. Table 6: Background-only fit results for the A- and B-type regions (top table) and C-type regions (bottom table) performed using 139  fb -1 of data. The quoted uncertainties on the fitted SM background include both the statistical and systematic uncertainties.

	SRA			SRA-L			SRA-M			SRA-H			SRB
Observed events	$17$			$12$			$3$			$2$			$3$
Fitted SM bkg events	$17.1$	$\pm$	$2.8$	$8.4$	$\pm$	$1.7$	$5.7$	$\pm$	$0.8$	$3.0$	$\pm$	$1.5$	$3.3$	$\pm$	$0.9$
$t\overline{t}$	$10.1$	$\pm$	$2.5$	$4.7$	$\pm$	$1.5$	$3.7$	$\pm$	$0.6$	$1.7$	$\pm$	$1.4$	$2.3$	$\pm$	$0.8$
$Z+\mathrm{jets}$	$2.6$	$\pm$	$0.4$	$1.3$	$\pm$	$0.2$	$0.9$	$\pm$	$0.2$	$0.4$	$\pm$	$0.1$	$0.3$	$\pm$	$0.1$
Single-top	$1.4$	$\pm$	$0.3$	$0.4$	$\pm$	$0.1$	$0.3$	$\pm$	$0.1$	$0.6$	$\pm$	$0.2$	$0.5$	$\pm$	$0.1$
$t\overline{t}$ +$W$/$Z$	$1.2$	$\pm$	$0.3$	$0.7$	$\pm$	$0.1$	$0.3$	$\pm$	$0.1$	$0.1$	$\pm$	$0.0$	$0.07$	$\pm$	$0.02$
$t\overline{t}$ +$h$	$1.1$	$\pm$	$0.2$	$0.7$	$\pm$	$0.1$	$0.3$	$\pm$	$0.1$	$0.1$	$\pm$	$0.0$	$0.13$	$\pm$	$0.02$
$W+\mathrm{jets}$	$0.4$	$\pm$	$0.1$	$0.2$	$\pm$	$0.1$	$0.1$	$\pm$	$0.0$		–		$0.02$	$\pm$	$0.01$
Diboson	$0.4$	$\pm$	$0.1$	$0.3$	$\pm$	$0.1$		–			–			–
$m({\tilde{b}}_1,{\tilde{\chi }}_2^0,{\tilde{\chi }}_1^0)=(1100,330,200)$ GeV	$13.7$	$\pm$	$0.3$	$0.7$	$\pm$	$0.1$	$6.3$	$\pm$	$0.2$	$6.6$	$\pm$	$0.2$	$0.3$	$\pm$	$0.1$
$m({\tilde{b}}_1,{\tilde{\chi }}_2^0,{\tilde{\chi }}_1^0)=(700,680,550)$ GeV	$1.3$	$\pm$	$0.6$	$0.2$	$\pm$	$0.1$	$0.5$	$\pm$	$0.4$	$0.6$	$\pm$	$0.4$	$7.4$	$\pm$	$1.2$
$m({\tilde{b}}_1,{\tilde{\chi }}_2^0,{\tilde{\chi }}_1^0)=(1200,1150,60)$ GeV	$8.7$	$\pm$	$0.2$	$1.4$	$\pm$	$0.1$	$3.4$	$\pm$	$0.1$	$3.8$	$\pm$	$0.1$	$0.6$	$\pm$	$0.1$

Table 7. Table 7: From left to right, observed 95% CL upper limits on the visible cross sections σ vis subscript 𝜎 vis \sigma_{\textrm{vis}} , the observed ( S obs 95 subscript superscript 𝑆 95 obs S^{95}_{\mathrm{obs}} ) and expected ( S exp 95 subscript superscript 𝑆 95 exp S^{95}_{\mathrm{exp}} ) 95% CL upper limits on the number of signal events with ± plus-or-minus \pm 1 σ 𝜎 \sigma excursions of the expectation, the CL of the background-only hypothesis, CL B subscript CL B \mathrm{CL}_{\mathrm{B}} , the discovery p 𝑝 p -value ( p 0 subscript 𝑝 0 p_{0} ), truncated at 0.5, and the associated significance (in parentheses).

Signal channel	${\sigma }_{\mathrm{vis}}$[fb]	$S_{\mathrm{obs}}^{95}$	$S_{\exp }^{95}$	${\mathrm{CL}}_B$	$p_0$ ($Z$)
SRA	$0.06$	$9.0$	${10.1}_{-3.1}^{+4.7}$	$0.38$	$0.50$ $\left(0.00\right)$
SRB	$0.04$	$4.9$	${5.1}_{-1.7}^{+2.8}$	$0.45$	$0.50$ $\left(0.00\right)$
SRC	$0.19$	$26.0$	${20.8}_{-5.5}^{+7.0}$	$0.80$	$0.17$ $\left(0.97\right)$

Equations

$$m_{\mathrm{eff}}=\sum _{i\le N_{\mathrm{jets}}}(p_{\mathrm{T}}^{\mathrm{jet}}{)}_i+E_{\text{T}}^{\text{miss}}.$$

$$\mathcal{S}=\sqrt{\frac{|{\mathbf{p}}_{\mathrm{T}}^{\mathrm{miss}}{|}^2}{{\sigma }_{\text{L}}^2(1-{\rho }_{\text{LT}}^2)}}.$$

Full text

\AtlasTitle

Search for bottom-squark pair production with the ATLAS detector in final states containing Higgs bosons, $b$-jets and missing transverse momentum \AtlasAbstract The result of a search for the pair production of the lightest supersymmetric partner of the bottom quark ($\tilde{b}_{1}$) using 139\,\mbox{fb{}^{-1}} of proton–proton data collected at $\sqrt{s}=13$ TeV by the ATLAS detector is reported. In the supersymmetric scenarios considered both of the bottom-squarks decay into a $b$-quark and the second-lightest neutralino, $\tilde{b}_{1}\rightarrow b+\tilde{\chi}^{0}_{2}$. Each $\tilde{\chi}^{0}_{2}$ is assumed to subsequently decay with 100% branching ratio into a Higgs boson ($h$) like the one in the Standard Model and the lightest neutralino: $\tilde{\chi}^{0}_{2}\rightarrow h+\tilde{\chi}^{0}_{1}$. The $\tilde{\chi}^{0}_{1}$ is assumed to be the lightest supersymmetric particle (LSP) and is stable. Two signal mass configurations are targeted: the first has a constant LSP mass of 60 GeV; and the second has a constant mass difference between the $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$ of 130 GeV. The final states considered contain no charged leptons, three or more $b$-jets, and large missing transverse momentum. No significant excess of events over the Standard Model background expectation is observed in any of the signal regions considered. Limits at the 95% confidence level are placed in the supersymmetric models considered, and bottom-squarks with mass up to 1.5 TeV are excluded.

\AtlasRefCodeSUSY-2018-031 \PreprintIdNumberCERN-EP-2019-142 \AtlasDate

\AtlasJournalRefJHEP 12 (2019) 060 \AtlasDOI10.1007/JHEP12(2019)060 \AtlasCoverSupportingNoteANA-SUSY-2018-31-INThttps://cds.cern.ch/record/2638536 \AtlasCoverCommentsDeadline9 July 2019 \AtlasCoverAnalysisTeamJohn Anders, Davide Costanzo, Monica D’Onofrio, Evangelos Kourlitis, Calum Macdonald, Hamish Teagle, Thomas Weston \AtlasCoverEdBoardMemberStan Lai (chair), Rui Wang, Stefano Passagio \AtlasCoverEgroupEditorsatlas-ana-susy-2018-31-analysis-team@cern.ch \AtlasCoverEgroupEdBoardatlas-ana-susy-2018-31-editorial-board@cern.ch

††journal: JHEP

1 Introduction

Supersymmetry (SUSY) [1, 2, 3, 4, 5, 6] provides an extension to the Standard Model (SM) that solves the hierarchy problem [7, 8, 9, 10] by introducing partners of the known bosons and fermions. In $R$-parity-conserving models [11], SUSY particles are produced in pairs and the lightest supersymmetric particle (LSP) is stable and provides a candidate for dark matter [12, 13]. The superpartners of the SM bosons (the wino, bino and higgsinos) mix to form the neutralinos ($\tilde{\chi}_{1,2,3,4}^{0}$) and charginos ($\tilde{\chi}_{1,2}^{\pm}$) physical states. For a large selection of models, the LSP is the lightest neutralino ($\tilde{\chi}^{0}_{1}$). Naturalness considerations suggest that the supersymmetric partners of the third-generation quarks are light [14, 15]. If this is assumed, the lightest bottom-squark ($\tilde{b}_{1}$) and lightest top-squark ($\tilde{t}_{1}$) mass eigenstates111The scalar partners of the left-handed and right-handed chiral components of the bottom quark ($\tilde{b}_{\mathrm{L},\mathrm{R}}$) or top quark ($\tilde{t}_{\mathrm{L},\mathrm{R}}$) mix to form mass eigenstates for which $\tilde{b}_{1}$and $\tilde{t}_{1}$are defined as the lighter of the two states. could be significantly lighter than the other squarks and the gluinos. As a consequence, $\tilde{b}_{1}$ and $\tilde{t}_{1}$ could be pair-produced with relatively large cross-sections at the Large Hadron Collider (LHC). Depending on the mass hierarchy considered, it is possible that the $\tilde{b}_{1}$ and $\tilde{t}_{1}$ could decay into final states with Higgs bosons, h, like the one in the SM, and this allows the Higgs boson to be used as a probe for new physics.

This article presents a search for the pair production of bottom squarks decaying into the LSP via a complex decay chain containing the second-lightest neutralino ($\tilde{\chi}^{0}_{2}$) and the Higgs boson: $\tilde{b}_{1}\rightarrow b+\tilde{\chi}^{0}_{2}$ and subsequently $\tilde{\chi}^{0}_{2}\rightarrow h+\tilde{\chi}^{0}_{1}$. Such a decay hierarchy is predicted in minimal supersymmetric extensions to the SM (MSSM) [16, 17], with $h$ assumed to be the lightest of the neutral bosons introduced in the MSSM. The bottom squark decaying through a next-to-lightest neutralino is one of the possible modes within the MSSM. Dedicated searches for direct decays into the lightest neutralino ($\tilde{b}_{1}\rightarrow b\tilde{\chi}^{0}_{1}$) or a chargino ($\tilde{b}_{1}\rightarrow t\tilde{\chi}_{1}^{\pm}$) have been reported by the ATLAS and CMS collaborations (see for example  [18, 19] and  [20, 21, 22]).

When the LSP is bino-like and the $\tilde{\chi}^{0}_{2}$ is a wino–higgsino mixture, the branching ratio ($\mathcal{B}$) of $\tilde{\chi}^{0}_{2}\rightarrow h+\tilde{\chi}^{0}_{1}$ is enhanced relative to the other possible $\tilde{\chi}^{0}_{2}$ decays. The Higgs boson mass is taken to be 125 GeV, and the decay into a pair of bottom quarks is assumed to be the same as in the SM ($\mathcal{B}$ = 58% [23, 24]), although it could be enhanced or reduced in the MSSM.

This search is interpreted within simplified model scenarios [25, 26] and Figure 1 illustrates the targeted model. In the first set of models, already considered by the ATLAS Collaboration using 8 TeV data [27], the mass of the $\tilde{\chi}^{0}_{1}$ is fixed at 60 GeV. The bottom-squark and $\tilde{\chi}^{0}_{2}$ masses vary in the ranges 250–1600 GeV and 200–1500 GeV, respectively. The assumption about the $\tilde{\chi}^{0}_{1}$ mass is motivated by dark-matter relic density measurements and might be favoured in Higgs-pole annihilation scenarios [28] where $m_{\tilde{\chi}^{0}_{1}}\simeq m_{h}/2$. The previous search performed by ATLAS using 8 TeV data excluded bottom-squark masses up to 750 GeV in this scenario [27].

The second set of SUSY models assumes a fixed mass difference between the $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$, sufficient to produce an on-shell Higgs boson. The mass difference, $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})$, is set to 130 GeV, whilst bottom-squark and $\tilde{\chi}^{0}_{1}$ masses vary in the ranges 400–1500 GeV and 1–800 GeV, respectively. A similar scenario is considered by the CMS Collaboration in Ref. [29], where the $h\rightarrow\gamma\gamma$ decay mode is exploited to exclude bottom-squark masses up to 530 GeV; no prior ATLAS searches have targeted these models.

The final states are characterised by a unique signature, which contains many jets, of which up to six can be identified as originating from the fragmentation of $b$-quarks (referred to as $b$-jets), missing transverse momentum ($\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$, the magnitude thereof referred to as $E_{\mathrm{T}}^{\mathrm{miss}}$), and no charged leptons (referred to as leptons). New selections and dedicated procedures aiming to maximise the efficiency of reconstructing the Higgs boson candidates decaying into a $b$-quark pair are employed in this article. Section 2 presents a brief overview of the ATLAS detector, with Section 3 describing the data and simulated samples used in the analysis. The event reconstruction methods are explained in Section 4. An overview of the analysis strategy is presented in Section 5, with the background estimation strategy discussed in Section 6. The systematic uncertainties considered in the analysis are described in Section 7. Section 8 presents the results and interpretation thereof, with the conclusions presented in Section 9.

2 ATLAS detector

The ATLAS detector [30] is a multipurpose particle physics detector with a forward–backward symmetric cylindrical geometry and nearly 4$\pi$ coverage in solid angle.222ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point in the centre of the detector. The positive $x$-axis is defined by the direction from the interaction point to the centre of the LHC ring, with the positive $y$-axis pointing upwards, while the beam direction defines the $z$-axis. Cylindrical coordinates $(r,\phi)$ are used in the transverse plane, $\phi$ being the azimuthal angle around the $z$-axis. The component of momentum in the transverse plane is denoted by $p_{\mathrm{T}}$. The pseudorapidity $\eta$ is defined in terms of the polar angle $\theta$ by $\eta=-\ln\tan(\theta/2)$. Rapidity is defined as $y$ = $0.5\ln[(E+p_{z})/(E-p_{z})]$ where $E$ denotes the energy, and $p_{z}$ is the component of the momentum along the beam direction. The separation of two objects in $\eta$–$\phi$ space is given by $\Delta R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$. The inner tracking detector consists of pixel and silicon microstrip detectors covering the pseudorapidity region $|\eta|<2.5$, surrounded by a transition radiation tracker which enhances electron identification in the region $|\eta|<2.0$. Between Run 1 and Run 2, a new inner pixel layer, the insertable B-layer [31, 32], was added at a mean sensor radius of 3.3 cm. The inner detector is surrounded by a thin superconducting solenoid providing an axial 2 T magnetic field and by a fine-granularity lead/liquid-argon (LAr) electromagnetic calorimeter covering $|\eta|<3.2$. A steel/scintillator-tile calorimeter provides hadronic coverage in the central pseudorapidity range ($|\eta|<1.7$). The endcap and forward regions ($1.5<|\eta|<4.9$) of the hadronic calorimeter are made of LAr active layers with either copper or tungsten as the absorber material. An extensive muon spectrometer with an air-core toroidal magnet system surrounds the calorimeters. Three layers of high-precision tracking chambers provide coverage in the range $|\eta|<2.7$, while dedicated fast chambers allow triggering in the region $|\eta|<2.4$. The ATLAS trigger system consists of a hardware-based level-1 trigger followed by a software-based high-level trigger [33].

3 Data and simulated event samples

The data analysed in this study correspond to a total of 139\,\mbox{fb{}^{-1}} of proton–proton ($pp$) collision data collected by the ATLAS detector with a centre-of-mass energy of 13 TeV and a 25 ns proton bunch crossing interval in the period between 2015 and 2018. All detector subsystems were required to be operational during data taking. The average number of interactions per bunch crossing (pile-up) increased from $\langle\mu\rangle$ $=20$ (2015–2016 dataset) to $\langle\mu\rangle$ $=37$ (2018 dataset), with a highest $\langle\mu\rangle$ $=38$ (2017 dataset). The uncertainty in the combined 2015–2018 integrated luminosity is 1.7 % [34], obtained using the LUCID-2 detector [35] for the primary luminosity measurements.

Events are required to pass an $E_{\mathrm{T}}^{\mathrm{miss}}$ trigger [36] which is fully efficient for events with reconstructed $E_{\mathrm{T}}^{\mathrm{miss}}>250\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$. Additional single-lepton triggers requiring electrons or muons are used to estimate the SM backgrounds, with an offline selection of $p_{\mathrm{T}}(\ell)>27$ GeV used to ensure the trigger is fully efficient ($\ell=e,\mu$).

Dedicated Monte Carlo (MC) simulated samples are used to model SM processes and estimate the expected signal yields. All samples were produced using the ATLAS simulation infrastructure [37] and GEANT4 [38], or a faster simulation based on a parameterisation of the calorimeter response and GEANT4 for the other detector systems [37].

The SUSY signal samples were generated with MadGraph5_aMC@NLO v2.6.2 [39] at leading order (LO) and interfaced to PYTHIA v8.230 [40] for the modelling of the parton showering (PS), hadronisation and the underlying event with the A14 [41] set of tuned parameters (tune). The matrix element (ME) calculation was performed at tree level and includes the emission of up to two additional partons. The ME–PS matching was done using the CKKW-L [42] prescription, with a matching scale set to one quarter of the bottom-squark mass. The NNPDF2.3 LO [43] parton distribution function (PDF) set was used. Signal cross-sections were calculated to approximate next-to-next-to-leading order in the strong coupling constant, adding the resummation of soft gluon emission at next-to-next-to-leading-logarithm (approximate NNLO+NNLL) [44, 45, 46, 47] accuracy. The nominal cross-section and its uncertainty were derived using the PDF4LHC15_mc PDF set, following the recommendations of Ref. [48]. For $\tilde{b}_{1}$ masses between 400 GeV and 1.5 TeV, the cross-sections range from 2.1 pb to 0.26 fb, with uncertainties from 6% to 17%.

The SM backgrounds considered in this analysis are: $t\bar{t}$ pair production; single-top-quark production; $Z\!+\!\mathrm{jets}$; $W\!+\!\mathrm{jets}$; $t\bar{t}$ production with an electroweak ($ttV$) or Higgs ($ttH$) boson; and diboson production. The samples were simulated using different MC generator programs depending on the process. Pair production of top quarks, $t\bar{t}$, was generated using POWHEG-BOX v2 [49, 50, 51, 52] interfaced with PYTHIA v8.230 and the A14 tune with the NNPDF2.3 LO PDF set for the ME calculations. The $h_{\mathrm{damp}}$ parameter in POWHEG-BOX, which controls the $p_{\mathrm{T}}$ of the first additional emission beyond the Born level and thus regulates the $p_{\mathrm{T}}$ of the recoil emission against the $t\bar{t}$ system, was set to 1.5 times the top-quark mass ($m_{t}\leavevmode\nobreak\ =\leavevmode\nobreak\ 172.5$ GeV) as a result of studies documented in Ref. [53]. The generation of single top quarks in the $Wt$-channel, $s$-channel and $t$-channel production modes was performed by POWHEG-BOX v2 [54, 50, 51, 52] similarly to the $t\bar{t}$ samples. For all processes involving top quarks, top-quark spin correlations were preserved. All events with at least one leptonically decaying $W$ boson were retained; fully hadronic $t\bar{t}$ and single-top events do not contain sufficient $E_{\mathrm{T}}^{\mathrm{miss}}$ to contribute significantly to the background. The production of $t\bar{t}$ pairs in association with electroweak vector bosons ($W,Z$) or Higgs bosons was modelled by samples generated at NLO using MadGraph5_aMC@NLO v2.2.3 and showered with PYTHIA v8.212. Events containing $W$ or $Z$ bosons with associated jets, including jets from the fragmentation of heavy-flavour quarks, were simulated using the SHERPA v2.2.1 [55] generator. Matrix elements were calculated for up to two additional partons at NLO and four partons at LO using the Comix [56] and OpenLoops [57] ME generators and were merged with the SHERPA PS [58] using the ME+PS@NLO prescription [59]. The NNPDF3.0 NNLO [43] PDF set was used in conjunction with a dedicated PS tune developed by the SHERPA authors. Diboson processes were also simulated using the SHERPA generator using the NNPDF3.0 NNLO PDF set. They were calculated for up to one ($ZZ$) or zero ($WW,WZ$) additional partons at NLO and up to three additional partons at LO. Other potential sources of backgrounds, such as the production of three or four top quarks or three gauge bosons, are found to be negligible. Finally, contributions from multijet background are estimated from data using a jet smearing procedure described in Ref. [60] and are found to be negligible in all regions.

All background processes are normalised to the best available theoretical calculation for their respective cross-sections. The NLO $t\bar{t}$ inclusive production cross-section is corrected to the theory prediction at NNLO in QCD including the resummation of NNLL soft-gluon terms calculated using Top++2.0 [61, 62, 63, 64, 65, 66, 67]. Samples of single-top events are normalised to the NLO cross-sections reported in Refs. [68, 69, 70].

For all samples, except those generated using SHERPA, the EvtGen v1.2.0 [71] program was used to simulate the properties of the bottom- and charm-hadron decays. All simulated events include a modelling of contributions from pile-up by overlaying minimum-bias $pp$ interactions from the same (in-time pile-up) and nearby (out-of-time pile-up) bunch crossings simulated in PYTHIA v8.186 and EvtGen v1.2.0 with the A3 [72] tune and the NNPDF2.3 LO set [43].

4 Event reconstruction

This search is based upon a selection of events with many $b$-jets, large missing transverse momentum and no charged leptons (electrons and muons) in the final state. All events are required to have a reconstructed primary vertex which is consistent with the beamspot envelope and consists of at least two associated tracks in the inner detector with $p_{\mathrm{T}}$ $>$ 500 MeV. If more than one vertex passing the above requirements is found, the one with the largest sum of the squares of transverse momenta of associated tracks [73] is chosen.

Jet candidates are reconstructed from three-dimensional clusters of energy in the calorimeter [74] with the anti-$k_{t}$ jet algorithm [75, 76] using a radius parameter of 0.4. The application of a jet energy scale (JES) correction derived from data and simulation [77] is used to calibrate the reconstructed jets. A set of quality criteria is applied to identify jets which arise from non-collision sources or detector noise [78] and any event which contains a jet failing to satisfy these criteria is removed. Additional jets that arise from pile-up interactions are rejected by applying additional track-based selections to jets with $p_{\mathrm{T}}<120$ GeV and $|\eta|<2.4$ [79], and the jet momentum is corrected by subtracting the expected average energy contribution from pile-up using the jet area method [80]. Jets are classified as either ‘baseline’ or ‘signal’; baseline jets are required to have $p_{\mathrm{T}}>20$ GeV and $|\eta|<4.8$ whilst signal jets are selected after resolving overlaps with electrons and muons, as described below, and must pass tighter requirements of $p_{\mathrm{T}}>30$ GeV and $|\eta|<2.8$.

Signal jets are identified as $b$-jets if they are within $|\eta|<2.5$ and are tagged by a multivariate algorithm which uses a selection of inputs including information about the impact parameters of inner-detector tracks, the presence of displaced secondary vertices and the reconstructed flight paths of $b$- and $c$-hadrons inside the jet [81]. The $b$-tagging algorithm used has an efficiency of 77%, determined in a sample of simulated $t\bar{t}$ events. It was chosen as part of the optimisation procedure and the corresponding misidentification rate is 20% for $c$-jets and 0.9% for light-flavour jets. To compensate for differences between data and MC simulation in the $b$-tagging efficiencies and mis-tag rates, correction factors are derived from data and applied to the samples of simulated events; details are found in Ref. [81].

Electron candidates are reconstructed from energy clusters in the electromagnetic calorimeter matched to a track in the inner detector and are required to satisfy a set of ‘loose’ quality criteria [82]. They are also required to lie within the fiducial volume $|\eta|<2.47$ and have $p_{\mathrm{T}}>4.5$ GeV. Muon candidates are reconstructed by matching tracks in the inner detector with tracks in the muon spectrometer. Muon candidates which have a transverse (longitudinal) impact parameter relative to the primary vertex larger than 0.2 mm (1 mm) are rejected to suppress muons from cosmic rays. Muon candidates are also required to satisfy ‘medium’ quality criteria [83] and have $|\eta|<2.5$ and $p_{\mathrm{T}}>4$ GeV. Electron (muon) candidates are matched to the primary vertex by requiring the transverse impact parameter ($d_{0}$) to satisfy $|d_{0}|/\sigma(d_{0})<$ 5 (3), and the longitudinal impact parameter ($z_{0}$) to satisfy $|z_{0}\sin\theta|<$ 0.5 mm. Lepton candidates remaining after resolving overlaps with baseline jets are called ‘baseline’ leptons. In the control regions where tighter lepton identification is required, ‘signal’ leptons are chosen from the baseline set with $p_{\mathrm{T}}$ $>$ 20 GeV and are required to be isolated from other activity in the detector using a criterion designed to accept at least 95% of leptons from $Z$ boson decays; details are found in Ref. [84]. In the dilepton control region where the single-lepton triggers are used, the leading lepton is required to have $p_{\mathrm{T}}$ $>$ 27 GeV; which ensures full efficiency of the single-lepton triggers. Signal electrons are further required to satisfy ‘tight’ quality criteria [82]. The MC events are corrected to account for differences in the lepton trigger, reconstruction and identification efficiencies between data and MC simulation.

Possible reconstruction ambiguities between baseline electrons, muons and jets are resolved by firstly removing electron candidates which share an inner detector track with a muon candidate. Jet candidates are then removed if they are within $\Delta R=\sqrt{(\Delta y)^{2}+(\Delta\phi)^{2}}=0.2$ of an electron candidate; next, electron candidates are discarded if they are within $\Delta R=0.4$ of a jet. Muons are discarded if they lie within $\Delta R=0.4$ of any remaining jet, except for the case where the number of tracks associated with the jet is less than three, where the muon is kept and the jet is discarded.

Identified $\tau$ leptons decaying hadronically are not considered but the following $\tau$-veto procedure is applied to reject events which contain $\tau$-like objects. Candidates ($\tau_{\mathrm{cand}}$) are identified as jets which have $|\eta|<2.5$ and less than five inner detector tracks of $p_{\mathrm{T}}>500$ MeV. If an event contains a tau candidate with a small azimuthal distance to the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$ ($\Delta\phi(E_{\mathrm{T}}^{\mathrm{miss}},\tau_{\mathrm{cand}})<\pi/3$), then the event is vetoed.

The missing transverse momentum $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$ is defined as the negative vector sum of the $p_{\mathrm{T}}$ of all selected and calibrated physics objects (electrons, muons, photons [82] and jets) in the event, with an extra term added to account for soft energy in the event which is not associated with any of the selected objects [85]. This soft term is calculated from inner-detector tracks with $p_{\mathrm{T}}$ above 500 MeV matched to the PV, thus ensuring it is robust against pile-up contamination [86, 87].

5 Analysis strategy

Three sets of non-orthogonal signal regions (SRs) are defined to target different mass hierarchies of the SUSY particles involved. These definitions exploit various discriminating observables and algorithms developed to explicitly reconstruct Higgs boson candidates in the decay chain. Events with charged leptons are vetoed in all SRs. Events with one or two charged leptons are used to define control regions (CRs) to aid in the estimation of the main SM backgrounds. Additionally, events with zero charged leptons are utilised to define validation regions (VRs) to ensure the background estimation method, described in Section 6, is robust. The optimisation procedure for the event selection aims to maximise the yield of bottom-squark pair production events while reducing SM background contributions. It is performed for the two simplified model scenarios introduced in Section 1. Since the $h\rightarrow bb$ decay mode is considered, the final state contains a large jet multiplicity, with many of these jets originating from $b$-quarks, and large $E_{\mathrm{T}}^{\mathrm{miss}}$ from the neutralinos.

The event selection criteria are defined on the basis of kinematic requirements for the objects described in the previous section and the event variables described below. For these definitions, signal jets are used and are ordered according to decreasing $p_{\mathrm{T}}$.

• •

$N_{\mathrm{jets}}$: the number of signal jets.

• •

$N_{b\textrm{-jets}}$: the number of $b$-jets.

• •

$\min\Delta\phi({\mathrm{jet}}_{1-4},\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}})$: the minimum azimuthal distance between the four highest-$p_{\mathrm{T}}$ jets and the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$. This is a powerful discriminating variable against multijet background events containing a large amount of $E_{\text{T}}^{\text{miss}}$ due to mismeasured jets. Typically, multijet background events exhibit low values of this variable and studies using data-driven multijet estimates indicate that a selection of $\min\Delta\phi({\mathrm{jet}}_{1-4},\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}})>0.4$ is sufficient to reduce the multijet background to a negligible level.

• •

$\Delta\phi(j_{1},\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}})$ : the azimuthal distance between the highest-$p_{\mathrm{T}}$ jet and the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$. This variable is used to select events where the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$ is expected to be recoiling against the leading jet.

• •

$m_{\mathrm{eff}}$: the effective mass [88] of an event is defined as the scalar sum of the $p_{\mathrm{T}}$ of all signal jets and the $E_{\text{T}}^{\text{miss}}$, i.e.:

[TABLE]

• •

$\mathcal{S}$: referred to as the “object-based $E_{\text{T}}^{\text{miss}}$-significance” [89] is defined as follows:

[TABLE]

The total momentum resolution of all jets and leptons, at a given $p_{\mathrm{T}}$ and $|\eta|$, is determined from parameterised Monte Carlo simulation which well reproduces the resolution measured in data. $\sigma_{\textup{L}}$ is the total momentum resolution after being rotated into the longitudinal (parallel to the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$) plane. The quantity $\rho_{\textup{LT}}$ is a correlation factor between the longitudinal and transverse momentum resolution (again with respect to the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$) of each jet or lepton. The significance $\mathcal{S}$ is used to discriminate events where the $E_{\mathrm{T}}^{\mathrm{miss}}$ arises from invisible particles in the final state from events where the $E_{\mathrm{T}}^{\mathrm{miss}}$ arises from poorly measured particles (and jets). Additionally, it is useful in discriminating between signal events with large $E_{\mathrm{T}}^{\mathrm{miss}}$ and $Z\!+\!\mathrm{jets}$ events with medium-to-low $E_{\mathrm{T}}^{\mathrm{miss}}$.

Additional selections on the $p_{\mathrm{T}}$ of the leading jet and of the leading $b$-jet are also applied as detailed in the following subsections. In all signal regions, events containing baseline leptons with $p_{\mathrm{T}}>10$ GeV are vetoed, as well as events containing $\tau$-lepton candidates that align with the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$ within $\Delta\phi=\pi/3$ . Only events with $E_{\mathrm{T}}^{\mathrm{miss}}$ $>250$ GeV are retained to ensure full efficiency of the trigger.

The event kinematics targeted by the three SRs are depicted in Figure 2. The first signal region is SRA, designed to target the ‘bulk’ region of both signal models, with moderate- to high-mass splitting between the $\tilde{b}_{1}$ and $\tilde{\chi}^{0}_{2}$. In these scenarios all of the $b$-jets, from both the bottom-squark and Higgs boson decays, are at a relatively high $p_{\mathrm{T}}$ and can be resolved in the detector. The $b$-jets from the Higgs boson can be isolated by removing the ones most likely from the bottom-squark decays and checking the angular separation between the remaining $b$-jets.

The second region, SRB, is designed to target the phase space of the $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=130$ GeV scenario with a small mass splitting between the $\tilde{b}_{1}$ and $\tilde{\chi}^{0}_{2}$, referred to as the “compressed” region. An initial-state radiation (ISR)-like selection is used where the small mass splitting between the bottom squark and neutralino leads to relatively soft $b$-jets from the bottom squark decay, which are difficult to reconstruct. In this scenario it is possible to reconstruct both Higgs bosons using angular separation methods. Finally, SRC is designed to target the “compressed” region of the $m(\tilde{\chi}^{0}_{1})=60$ GeV signal scenario, where the mass splitting between the $\tilde{b}_{1}$ and $\tilde{\chi}^{0}_{2}$ is small. The $b$-jets from the bottom squark decay are very soft and as such a lower $b$-jet multiplicity is used in this region, when compared to the A- and B-type selections. Additionally, the visible system ($b$-jets from the bottom squark decay and Higgs boson decay) is produced back-to-back with the reconstructed $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$.

5.1 The SRA selections

To exploit the kinematic properties of the signal over a large range of $\tilde{b}_{1}$, $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$ masses, incremental thresholds are imposed on the main discriminating variable, $m_{\mathrm{eff}}$, resulting in three mutually exclusive regions, 1.0 $<$ $m_{\mathrm{eff}}$ $<$ 1.5 TeV, 1.5 $<$ $m_{\mathrm{eff}}$ $<$ 2.0 TeV and $m_{\mathrm{eff}}>$ 2.0 TeV. These are labelled as SRA-L, -M and -H, respectively, to maximise coverage across the $\tilde{b}_{1}$ mass range. The selection criteria for the three SRAs are summarised in Table 1.

At least four $b$-tagged jets are required. To discriminate against multijet background, events where the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$ is aligned with a jet in the transverse plane are rejected by requiring $\min\Delta\phi({\mathrm{jet}}_{1-4},\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}})>0.4$. As a large $E_{\mathrm{T}}^{\mathrm{miss}}$ is expected from the neutralinos which escape the detector, a selection of $E_{\mathrm{T}}^{\mathrm{miss}}>350$ GeV is used. Additionally, the leading $b$-jet ($b_{1}$) is expected to have a large $p_{\mathrm{T}}$, hence a selection of $p_{\mathrm{T}}(b_{1})>200$ GeV is employed. At least one of the two Higgs boson candidates in the event is identified using a reconstruction algorithm referred to as max-min, which is a two-step procedure to remove the high-$p_{\mathrm{T}}$ $b$-jets from the bottom squark decay and then use the remaining $b$-jets to reconstruct a Higgs boson in the decay chain. The procedure is implemented as follows: first, pairs of $b$-jets are formed by iterating through all of the $b$-jets in the event, and the pair with the largest separation in $\Delta R$ is designated as arising from the bottom-squark decay and removed from the subsequent step; second, the pair with the smallest $\Delta R$ is identified as a possible Higgs boson candidate and its invariant mass calculated. The following $\Delta R$ and mass quantities are defined:

• •

$\Delta R_{\mathrm{max}}(b,b)$: the distance in $\eta$–$\phi$ between the two $b$-jets with the maximal angular separation which are most likely to originate from the initial decay of the $\tilde{b}_{1}$;

• •

$\Delta R_{\textrm{max-min}}(b,b)$: the distance in $\eta$–$\phi$ between the two $b$-jets with the minimum angular separation which are most likely to originate from the same Higgs boson decay, selected out of the remaining $b$-jets;

• •

$m(h_{\mathrm{cand}})$: the invariant mass of the $b$-jet pair identified as a Higgs candidate by the max-min algorithm. A lower bound on $m(h_{\mathrm{cand}})$ is used; in the majority of events the distribution peaks around the Higgs boson mass, but in scenarios where the incorrect combination of $b$-jets is chosen the signal can extend to higher masses.

When applied to signal, the max-min algorithm correctly selects a $h\rightarrow bb$ pairing in 20%–40% of cases for a single Higgs boson decay, depending upon the model. For a signal model corresponding to $m(\tilde{b}_{1},\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=(1100,330,200)$ GeV, about 3% of the simulated signal events are retained by the SRA selections.

5.2 The SRB selections

The SRB region targets small mass-splitting between the $\tilde{b}_{1}$and $\tilde{\chi}^{0}_{2}$ (of order 5–20 GeV), in the case of the $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=130$ GeV scenarios. The presence of an ISR jet boosting the bottom squarks, and consequently their decay products, is exploited. To efficiently suppress SM background contributions, events are selected where the highest-$p_{\mathrm{T}}$ jet is not $b$-tagged and has $p_{\mathrm{T}}$ $>$ 350 GeV; this jet is presumed to arise from ISR in the scenario under consideration. Additional selections of $E_{\mathrm{T}}^{\mathrm{miss}}$ $>350$ GeV and $\Delta\phi(j_{1},E_{\text{T}}^{\text{miss}})>2.8$ are applied. An $m_{\mathrm{eff}}$ selection of $>$ 1 TeV is also applied. The soft $p_{\mathrm{T}}$ spectrum predicted for $b$-jets from $\tilde{b}_{1}$decays can cause the $b$-jets to be difficult to reconstruct, hence a different algorithm, aiming to reconstruct both Higgs boson candidates, is employed.

Differently from the scenarios targeted by SRA, pairs of $b$-jets with the largest $\Delta R$ are found to be more likely to arise from the decay of the same Higgs boson candidate. Two pairs at a time are identified following an iterative procedure, such that at first the pair of $b$-jets leading to the highest $\Delta R$, $\Delta R_{bb1}$, is defined, followed by the second highest $\Delta R$, $\Delta R_{bb2}$, built considering only the remaining $b$-jets. The average mass of the two candidates $m(h_{\mathrm{cand1}},h_{\mathrm{cand2}})_{\mathrm{avg}}$ is calculated and a requirement is placed on the average mass, corresponding to a window around the Higgs boson mass: [75, 175] GeV. For a signal model corresponding to $m(\tilde{b}_{1},\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=(700,680,550)$ GeV, about 0.1% of the simulated signal events are retained by the SRB selections. The efficiency of correctly selecting the $b$-jets using this algorithm is in the range 15%–30%. The SRB requirements are listed in Table 2.

5.3 The SRC selections

When considering the scenario with a constant $\tilde{\chi}^{0}_{1}$ mass of 60 GeV, the $\Delta R$-based Higgs boson reconstruction algorithms are ineffective in the compressed region of phase space with a small mass splitting between the $\tilde{b}_{1}$ and $\tilde{\chi}^{0}_{2}$. In the inclusive SRC, the main discriminating quantity is $\mathcal{S}$; a selection of $\mathcal{S}>22$ is employed. Events are also required to have at least three $b$-jets. Four non-overlapping regions (SRC22, SRC24, SRC26 and SRC28) are defined as subsets of the inclusive SRC region, with incremental thresholds placed on $\mathcal{S}$ as detailed in Table 3, to ensure full coverage of the target models as a function of bottom-squark and neutralino mass. For a signal model corresponding to $m(\tilde{b}_{1},\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})=(1200,1150,60)$ GeV, about 11% of the simulated signal events are retained by the SRC selections. The $\mathcal{S}$ variable is effective in rejecting the SM background arising from associated production of a $Z$ boson decaying into neutrinos and $b$-jets.

6 Background estimation

There are two main SM backgrounds which are expected to contribute to the yields for the SRs introduced in the previous section. For SRAs and SRB, the main background is top-quark production which, according to MC estimates, contributes between 70% and 85% of the total background, depending upon the region considered, and is dominated by top-quark pairs produced in association with two $b$-quarks arising from gluon splitting. In the SRCs, the main backgrounds arise from $Z\!+\!\mathrm{jets}$ (up to 50% of the total) and from top-quark-related processes (up to 20% of the total).

The main SM backgrounds in each SR are determined separately with a profile likelihood fit to the event yields in the associated CRs [90]. This is commonly referred to as a background-only fit which constrains and adjusts the normalisation of the background processes. The background-only fit uses the observed event yield and the expected number of MC events in the associated CRs, which are described by Poisson statistics, as a constraint to adjust the normalisation of the background processes assuming that no signal is present.

The normalisation factor is referred to as the $\mu$ factor. The CRs are designed to be enriched in specific background contributions relevant to the analysis, whilst minimising the potential signal contamination, and they are orthogonal to the SRs.

When performing the fit for SRA, a multi-bin approach is used, with a single CR divided into three bins of $m_{\mathrm{eff}}$. Such an approach allows the calculation and use of a single normalisation parameter (applied to the main $t\bar{t}$ background across all bins of $m_{\mathrm{eff}}$), and additionally enables the fit to take into account the modelling of the $m_{\mathrm{eff}}$ variable.

The systematic uncertainties, described in Section 7, are included in the fit as nuisance parameters. They are constrained by Gaussian distributions with widths corresponding to the sizes of the uncertainties and are treated as correlated, when appropriate, between the various regions. The product of the various probability density functions and the Gaussian distributions forms the likelihood function, which the fit maximises by adjusting the background normalisation and the nuisance parameters. This approach reduces the influence of systematic uncertainties on the backgrounds with dedicated CRs, as these are absorbed by the normalisation parameter.

Finally, the reliability of the MC extrapolation of the SM background estimates outside of the CRs is evaluated in dedicated VRs, orthogonal to CRs and SRs.

The fit strategies for the A- and B-type regions are very similar and are represented schematically in Figure 3(a). They rely on CRs with a single-lepton requirement, as the $t\bar{t}$ background in the SR is dominated by semileptonic $t\bar{t}$ decays where the lepton is not identified. The main background in both regions is $t\bar{t}$ pair production in association with heavy-flavour jets. The fit strategy for the C-type regions is presented in Figure 3(b). The strategy is different because the main background in these regions is $Z$+jets, closely followed by the top-quark backgrounds. In order to define CRs enhanced in $t\bar{t}$ and $Z$+jets, additional variables are used:

• •

$m_{\mathrm{T}}$ : the event transverse mass $m_{\mathrm{T}}$ is defined as $m_{\mathrm{T}}=\sqrt{2p_{\mathrm{T}}(\ell)E_{\mathrm{T}}^{\mathrm{miss}}(1-\cos(\Delta\phi))}$, where $\Delta\phi$ is the difference in azimuthal angle between the lepton and the $\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}}$. This is used in the one-lepton CRs to reject multi-jet events which can be misidentified as containing a prompt lepton.

• •

$m_{\ell\ell}$: the invariant mass of the two leptons in the event. Since the two-lepton CR is used to constrain the $Z$+jets background, the $m_{\ell\ell}$ variable is required to be within the $Z$-mass window: [86, 106] GeV (used exclusively in the two-lepton CR).

• •

$\tilde{E}_{\mathrm{T}}^{\mathrm{miss}}$: the ‘lepton corrected’ $E_{\text{T}}^{\text{miss}}$. For the two-lepton CR the transverse momentum vectors of the leptons are subtracted from the $E_{\text{T}}^{\text{miss}}$ calculation in order to mimic the neutrinos from $Z\rightarrow\nu\nu$ decays (used exclusively in the two-lepton CR).

When designing the CRs and VRs, the potential signal contamination is checked in each region to ensure that the contribution from the signal process being targeted is small in the regions. The signal contamination in the CRs and VRs is found to be negligible, at the level of $<$ 1% of the total SM expectation, depending upon the signal mass hierarchy of the models considered in this search.

6.1 A-type CR and VR definitions

A single, $t\bar{t}$-dominated CR (CRA1$\ell$) is defined for the A-type regions and is split into the same three identical $m_{\mathrm{eff}}$ selections as the SRAs. The CR is defined similarly to the SR selection (as documented in Table 1); however, exactly one signal lepton (either $e$ or $\mu$) with $p_{\mathrm{T}}$ $>$ 20 GeV is required in the final state. Furthermore, the selections used to isolate the Higgs boson in the SRAs, namely the $\Delta R_{\mathrm{max}}(b,b)$, $\Delta R_{\textrm{max-min}}(b,b)$ and $m(h_{\mathrm{cand}})$ selections, are not applied in order to increase the number of events in the CR. The leading $b$-jet $p_{\mathrm{T}}$ selection is lowered to $>100$ GeV to further increase the number of events in the region, and a selection on the transverse mass of $m_{\mathrm{T}}>20$ GeV is applied to suppress misidentified leptons. Such selections result in pure $t\bar{t}$ CRs, with $t\bar{t}$ contributing more than 80% of the total SM contribution in each of the CRs. The fraction of top-quark pairs produced in association with $b$-quarks is equivalent between CRs and SRs, and accounts for about 70% of the total $t\bar{t}$ background. Figure 4(a) presents the distribution of $m(h_{\mathrm{cand}})$ in CRA1$\ell$, and shows that this variable is well modelled.

A zero-lepton validation region (VRA0$\ell$) is also defined, and split according to the same $m_{\mathrm{eff}}$ thresholds as the SRAs and CRAs. This VR is used to validate the modelling of the $t\bar{t}$ background when extrapolating from the one-lepton CRs to zero-lepton regions. The selections are based upon the SR selections but the VRs are orthogonal due to the $b$-jet multiplicity selection, which requires exactly three $b$-jets. Additionally, the $\Delta R_{\mathrm{max}}(b,b)$, $\Delta R_{\textrm{max-min}}(b,b)$ and $m(h_{\mathrm{cand}})$ selections are not applied in this region. A selection of $\mathcal{S}$ $<$ 22 is applied to ensure this region is orthogonal to the SRC regions.

6.2 B-type CR and VR definitions

For the B-type $t\bar{t}$ CR (CRB1$\ell$), a similar method of using a one-lepton region enriched in $t\bar{t}$ is implemented. The SR selections (as documented in Table 2) are applied, and additionally exactly one signal lepton with $p_{\mathrm{T}}$ $>$ 20 GeV is required. The $m(h_{\mathrm{cand1}},h_{\mathrm{cand2}})_{\mathrm{avg}}$ selection is dropped to increase the number of events in the region, and the $|\Delta\phi(j_{1},E_{\text{T}}^{\text{miss}})|$ selection is loosened to $>2.2$. Similarly to the A-type CR, a selection of $m_{\mathrm{T}}>20$ GeV is applied to suppress misidentified leptons. These selections result in a pure CR with 80% of the total expected SM background consisting of $t\bar{t}$. Figure 4(b) presents the $m(h_{\mathrm{cand1}},h_{\mathrm{cand2}})_{\mathrm{avg}}$ distribution in this region; it is shown to be well modelled.

The associated VR (VRB0$\ell$) is defined in a similar manner to the A-type VR, with selections similar to those of the SRB region, but an exclusive $b$-jet multiplicity selection of exactly three $b$-jets. Additionally, the selections used to reconstruct the Higgs bosons in the event are dropped to enhance the number of events in the region. A selection of $\mathcal{S}$ $<$ 22 is also applied to ensure this region is orthogonal to the C-type SRs.

6.3 C-type CR and VR definitions

Two CRs are defined for the C-type SRs, one to constrain the $Z$+jets background (CRC2$\ell$) and one to constrain the backgrounds associated with top quarks, $t\bar{t}$ and single top (CRC1$\ell$). A single normalisation parameter is used to constrain both the $t\bar{t}$ and single-top backgrounds, while the $Z\!+\!\mathrm{jets}$ background is constrained with an additional normalisation parameter. These CRs are based upon the SR shown in Table 3, but are orthogonal due to the different lepton multiplicities required.

The CRC2$\ell$ requires two same-flavour (SF) opposite-sign (OS) leptons, with invariant mass in the $Z$-mass window. The leading two leptons are required to have $p_{\mathrm{T}}$ $>$ 27 GeV and $p_{\mathrm{T}}$ $>$ 20 GeV respectively. To imitate the $E_{\text{T}}^{\text{miss}}$ selection in the SR, a selection of $\tilde{E}_{\mathrm{T}}^{\mathrm{miss}}>250$ GeV is utilised. For this region the selections on $\mathcal{S}$ are dropped to enhance the number of events in the region. Figure 4(c) shows the $\tilde{E}_{\mathrm{T}}^{\mathrm{miss}}$ distribution in this region. The CRC1$\ell$ region used to constrain the top-quark-related backgrounds requires one signal lepton with $p_{\mathrm{T}}$ $>$ 20 GeV. A selection of $\mathcal{S}$ $>$ 17 is applied. Similarly to the A- and B-type CRs, a selection of $m_{\mathrm{T}}$ $>$ 20 GeV is applied to remove the multi-jet contribution with fake or non-prompt leptons. Figure 4(d) presents the $\mathcal{S}$ distribution in this region.

Two zero-lepton VRs are defined to validate the extrapolation from CR to SR based on the SR selections. A VR with zero leptons and two $b$-jets (VRC0$\ell$-Z) with $\mathcal{S}$ $\in[20,22]$ and $E_{\text{T}}^{\text{miss}}$ $\in[250,600]$ GeV ensures a region orthogonal to the SR, but with a large contribution from the $Z$+jets process. A second VR is used to validate the modelling of the $t\bar{t}$ and single-top backgrounds (VRC0$\ell$-T); a selection of zero leptons, $\mathcal{S}$ $\in[15,22]$ and an inverted selection on the $\min\Delta\phi({\mathrm{jet}}_{1-4},\mathbf{p}_{\mathrm{T}}^{\mathrm{miss}})$ $\in[0.2,0.4]$ is applied to ensure orthogonality.

6.4 Summary of CR and VR results

A full overview of the control and validation regions used in the analysis can be found in Table 4. The control region pre-fit yields and fitted normalisation factors $\mu_{\mathrm{bkg}}$ for the A-, B- and C-type regions are presented in Figure 5(a). All $\mu$ values are consistent with unity, within 2$\sigma$ of the normalisation uncertainty, suggesting the modelling of the key SM background processes is already good before performing the fit. Figure 5(b) presents the observed yields, post-fit background estimates and significance [91] for the A-, B- and C-type validation regions. The background-only fit estimates are in good agreement with the data in these regions, and the post-fit expectation is within 1$\sigma$ of the central value for all regions.

7 Systematic uncertainties

Several sources of experimental and theoretical systematic uncertainty on the signal and background estimates are considered in this analysis. Their impact is reduced by fitting the event yields and normalising the dominant backgrounds in the CRs defined with kinematic selections resembling those of the corresponding SRs (see Section 6). Uncertainties due to the numbers of events in the CRs are also introduced in the fit for each region. The magnitude of the contributions arising from detector, theoretical modelling and statistical uncertainties are summarized in Table 5.

Dominant detector-related systematic uncertainties arise from the $b$-tagging efficiency and mis-tagging rates, and from the jet energy scale and resolution. In SRA and SRB, the contributions of these uncertainties are almost equivalent. In SRC, the $b$-tagging uncertainty is dominant. The systematic uncertainty on the $b$-tagging efficiency ranges from 4.5% for $b$-jets with $p_{\mathrm{T}}$ $\in[35,40]$ GeV up to 7.5% for $b$-jets with high $p_{\mathrm{T}}$ ($>100$ GeV). The $b$-tagging uncertainty is estimated by varying the $\eta$-, $p_{\mathrm{T}}$- and flavour-dependent scale factors applied to each jet in the simulation within a range that reflects the systematic uncertainty in the measured tagging efficiency and mis-tag rates in data  [81]. The uncertainties in the jet energy scale and resolution are based on their respective measurements in data [77, 92].

The uncertainties associated with lepton reconstruction and energy measurements have a negligible impact on the final results; however, the lepton, photon and jet-related uncertainties are propagated to the calculation of the $E_{\mathrm{T}}^{\mathrm{miss}}$, and additional uncertainties due to the energy scale and resolution of the soft term are included in the $E_{\mathrm{T}}^{\mathrm{miss}}$.

The systematic uncertainties related to the modelling of the energy of jets and leptons in the simulation are propagated to $\mathcal{S}$. No additional uncertainty on the energy resolution is applied, as the resolutions are taken to be the maximum of the parameterised data and simulation resolutions when performing the calculation for both data and MC simulation.

Uncertainties in the modelling of the SM background processes from MC simulation and their theoretical cross-section uncertainties are also taken into account. The dominant uncertainties in SRA and SRB arise from theoretical and modelling uncertainties of the $t\bar{t}$ background. They are computed as the difference between the predictions from nominal samples and those from additional samples differing in hard-scattering generator and parameter settings, or by using internal weights assigned to the events depending on the choice of renormalisation and factorisation scales, initial- and final-state radiation parameters, and PDF sets. The impact of the PS and hadronisation model is evaluated by comparing the nominal generator with a POWHEG sample interfaced to HERWIG 7 [93, 94], using the H7UE set of tuned parameters [94]. To assess the uncertainty due to the choice of hard-scattering generator and matching scheme, an alternative generator setup using aMC@NLO+PYTHIA8 is employed. It uses the shower starting scale, $\mu_{\mathrm{q}}=H_{\mathrm{T}}/2$, where $H_{\mathrm{T}}$ is defined here as the scalar sum of the $p_{\text{T}}$ of all outgoing partons.

The dominant uncertainties in SRC arise from the MC modelling of the $Z$+jets process, followed by the $t\bar{t}$ and single-top modelling. The $Z\!+\!\mathrm{jets}$ (as well as $W\!+\!\mathrm{jets}$) modelling uncertainties are estimated by considering different merging (CKKW-L) and resummation scales using alternative samples, PDF variations from the NNPDF3.0 NNLO replicas [55], and variations of factorisation and renormalisation scales in the ME. The latter have been evaluated using 7 point-variations, changing the renormalisation and factorisation scales up and down by factors 0.5 and 2, such that when one scale is up the other is down, and vice-versa.

For the SUSY signal processes, both the experimental and theoretical uncertainties in the expected signal yield are considered. Experimental uncertainties are found to be 6–36% across the mass plane with fixed LSP mass for A-type SRs, and 4–40% for C-type SRs. For models where $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})$ = 130 GeV is assumed, scenarios where SRB is relevant have uncertainties of 11–37%.

In all SRs, the dominant uncertainty on the signal yields is found to be from the $b$-tagging efficiency.

Theoretical uncertainties in the approximate NNLO+NNLL cross-section are calculated for each SUSY signal scenario, and are dominated by the uncertainties in the renormalisation and factorisation scales, followed by the uncertainty in the PDFs. These are 7–17% for bottom-squark masses in the range between 400 GeV and 1500 GeV. Additional uncertainties in the acceptance and efficiency due to the modelling of ISR and CKKW scale variations in SUSY signal MC samples are also taken into account, and contribute up to $\sim$10%.

8 Results and interpretation

The event yields for all SRs are reported in Table 6. The SM background expectations resulting from background-only fits are also reported showing statistical plus systematic uncertainties. The largest background contribution in A-type and B-type SRs arises from $t\bar{t}$  production, whilst the contribution from $Z\rightarrow\nu\bar{\nu}$ production in association with $b$-quarks is largest in the C-type SRs, with sub-dominant contributions from the $t\bar{t}$ and single-top processes. Other background sources are $t\bar{t}$ +$W/Z$, $t\bar{t}$ +$h$, diboson and $W$+jets production. The results are also summarised in Figure 6, where the significances for each of the SRs are also presented. No significant deviations are observed between expected and observed yields in all signal regions considered.

Figure 7 shows comparisons between the observed data and the post-fit SM predictions for some relevant kinematic distributions for the inclusive SRA, SRB and SRC selections before selection requirements are applied on the quantity shown. The expected distributions for scenarios with different bottom squark, $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$ masses (depending on the SR considered) are shown for illustrative purposes.

The CLs technique [95] is used to place 95% Confidence Level (CL) upper limits on event yields from physics beyond the SM (BSM) for each signal region. The profile-likelihood-ratio test statistic is used to exclude the signal-plus-background hypothesis for specific signal models. When normalised to the integrated luminosity of the data sample, results can be interpreted as corresponding upper limits on the visible cross-section, $\sigma_{\textrm{vis}}$, defined as the product of the BSM production cross-section, the acceptance and the selection efficiency of a BSM signal. When calculating the model-independent upper limits of the A- and C-type regions, only the inclusive SR selection is used. Table 7 summarises the observed ($S^{95}_{\mathrm{obs}}$) and expected ($S^{95}_{\mathrm{exp}}$) 95% CL upper limits on the number of BSM events and on $\sigma_{\textrm{vis}}$ for all SRs. The $p_{0}$-values, which represent the probability of the SM background to fluctuate to the observed number of events or higher, are also provided and are capped at $p_{0}=0.5$; the associated significance is provided in parentheses.

Model-dependent exclusion limits are obtained assuming the two types of SUSY particle mass hierarchies described in Section 1. The lightest bottom squark decays exclusively via $\tilde{b}_{1}\rightarrow b\tilde{\chi}^{0}_{2}$ with subsequent decay $\tilde{\chi}^{0}_{2}\rightarrow h\tilde{\chi}^{0}_{1}$. The expected limits from the SRs are compared for each set of scenarios and the observed limits are obtained by choosing the SR with the best expected sensitivity for each SUSY model. The fit procedure takes into account correlations in the yield predictions between control and signal regions due to common background normalisation parameters and systematic uncertainties. The experimental systematic uncertainties in the signal are taken into account for the calculation and are assumed to be fully correlated with those in the SM background.

Figures 8(a) and 8(b) show the observed (solid line) and expected (dashed line) exclusion contours at 95% CL in the $\tilde{b}_{1}$–$\tilde{\chi}^{0}_{2}$ mass planes for the two types of SUSY scenarios considered. For the scenarios where the mass of the neutralino is assumed to be 60 GeV, the sensitivity to models with the largest mass difference between the $\tilde{b}_{1}$ and the $\tilde{\chi}^{0}_{2}$ is achieved with the combination of the A-type SRs. Sensitivity to scenarios with small mass differences is obtained with the dedicated C-type SRs. For scenarios with $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})$ = 130 GeV, the sensitivity of the A-type SRs is complemented by the B-type SR in the case of small mass difference between the $\tilde{b}_{1}$ and the $\tilde{\chi}^{0}_{2}$.

Bottom-squark masses up to 1.5 TeV are excluded for models with fixed $m_{\tilde{\chi}^{0}_{1}}=60$ GeV and $\tilde{\chi}^{0}_{2}$ masses $\in[0.5,1.1]$ TeV. In case of $\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})$ = 130 GeV, bottom-squark masses up to 1.3 TeV are excluded for $\tilde{\chi}^{0}_{2}$ masses up to 750 GeV. The losses in sensitivity for models where $\tilde{\chi}^{0}_{2}$ masses are below 190 GeV are due to the stringent requirements on $E_{\mathrm{T}}^{\mathrm{miss}}$.

The results constitute a large improvement upon previous Run-1 searches and significantly strengthen the constraints on bottom squark masses; they are also complementary to other searches where bottom squarks are assumed to decay directly to a bottom-quark and a neutralino or to a top-quark and a chargino [96].

9 Conclusion

The result of a search for pair production of bottom squarks is reported. The analysis uses $139$ fb*-1* of $pp$ collisions at $\sqrt{s}=13$ TeV collected by the ATLAS experiment at the LHC between 2015 and 2018. $R$-parity-conserving SUSY scenarios where bottom squarks decay into a $b$-quark and the second-lightest neutralino, $\tilde{b}_{1}\rightarrow b+\tilde{\chi}^{0}_{2}$, with $\tilde{\chi}^{0}_{2}$ subsequently decaying into a Higgs boson like the one in the SM and the lightest neutralino, are considered. The search investigates final states containing large missing transverse momentum and three or more $b$-jets. No significant excess of events above the expected Standard Model background is found and exclusion limits at the 95% confidence level are placed on the visible cross-section and on the mass of the bottom squark for various assumptions about the mass hierarchy of the $\tilde{b}_{1}$, $\tilde{\chi}^{0}_{2}$ and $\tilde{\chi}^{0}_{1}$. Bottom-squark masses up to 1.5 (1.3) TeV are excluded for $\tilde{\chi}^{0}_{2}$ masses up to 1100 (750) GeV in models with fixed $m({\tilde{\chi}^{0}_{1}})=60$ GeV ($\Delta m(\tilde{\chi}^{0}_{2},\tilde{\chi}^{0}_{1})$ = 130 GeV). As the first search for such scenarios carried out by ATLAS in Run 2, these results are a significant improvement upon the previous Run-1 result, considerably tightening the constraints on bottom-squark production.

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