1 Introduction

Supersymmetry (SUSY) [16] provides an extension of the Standard Model (SM) that solves the hierarchy problem [710] by introducing supersymmetric partners of the known bosons and fermions. In the framework of the R-parity-conserving minimal supersymmetric extension of the SM (MSSM) [1113], SUSY particles are produced in pairs and the lightest supersymmetric particle (LSP) is stable, providing a possible candidate for dark matter [14, 15]. In a large variety of models the LSP is the lightest neutralino (\(\tilde{\chi }^{0}_{1}\)). Naturalness considerations [16, 17] suggest that the supersymmetric partners of the third-generation SM quarks are the lightest coloured supersymmetric particles. This may lead to the lightest bottom squark (\(\tilde{b}^{}_{1} \)) and top squark (\(\tilde{t}^{}_{1}\)) mass eigenstatesFootnote 1 being 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).

This paper presents a search for the pair production of bottom squarks decaying exclusively as \(\tilde{b}^{}_{1} \rightarrow b\tilde{\chi }^{0}_{1}\) using 3.2 fb\(^{-1}\)  of proton–proton (pp) collisions at \(\sqrt{s}=13\) TeV collected by the ATLAS experiment at the LHC in 2015. A SUSY particle mass hierarchy with this exclusive decay has been predicted by various phenomenological MSSM models [18]. Searches with the \(\sqrt{s}=8\) TeV LHC Run-1 dataset at ATLAS and CMS have set limits on \(\tilde{b}^{}_{1} \) masses in such scenarios. For \(\tilde{\chi }^{0}_{1}\) masses around 100 GeV, exclusion limits at 95 % confidence level (CL) up to 620 and 680 GeV have been reported by the ATLAS [19] and CMS [20] collaborations, respectively. The searches were performed in events characterized by the presence of large missing transverse momentum and two jets identified as containing b-hadrons (\(b\)-jets). In Run 2 of the LHC the production cross-section rises due to the increase of the centre-of-mass energy of the pp collisions. For instance, for a \(\tilde{b}^{}_{1} \) with a mass of 800 GeV, the production cross-section increases by almost a factor of ten going from \(\sqrt{s}=8\) TeV to \(\sqrt{s}=13\) TeV. In addition, the sensitivity of the analysis benefits from the improved algorithms adopted to identify b-jets and use of information from the newly installed pixel layer in the Run-2 ATLAS detector. The search strategy at 13 TeV closely follows the previous ATLAS studies, with signal regions defined to provide sensitivity to the kinematic topologies arising from different mass splittings between the bottom squark and the neutralino.

2 ATLAS detector

The ATLAS detector [21] is a multi-purpose particle physics detector with a forward-backward symmetric cylindrical geometry and nearly 4\(\pi \) coverage in solid angle.Footnote 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 (IBL) [22], was inserted 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 toroid 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 [23].

3 Data and simulated event samples

The data used in this analysis were collected by the ATLAS detector in pp collisions at the LHC with a centre-of-mass energy of 13 TeV and 25 ns proton bunch crossing interval during 2015. After applying beam-, data- and detector-quality criteria, the available dataset corresponds to an integrated luminosity of 3.2 fb\(^{-1}\) with an uncertainty of ±5 %. The uncertainty is derived from a calibration of the luminosity scale using a pair of xy beam-separation scans performed in August 2015 and following a methodology similar to that detailed in Ref. [24]. In this dataset, each event includes an average of approximately 14 additional inelastic pp collisions in the same bunch crossing (in-time pile-up). The events used in this search were selected using a trigger logic that accepts events with an uncorrected missing transverse momentum above 70 GeV, calculated using a sum over calorimeter cells. Additional events selected with lepton- and photon-based triggers are employed for control regions defined to estimate SM background contributions. All the selections employed in this paper use a highly efficient trigger selection. The triggers used are in a plateau region and the systematic uncertainties related to the trigger simulation are found to be negligible.

Monte Carlo (MC) simulated event samples are used to model the expected signal and to aid in the description and estimation of SM background processes. The response of the detector is simulated [25] either fully by a software program based on GEANT4  [26] or by a faster simulation based on a parameterization [27] for the calorimeter response and GEANT4 for the other detector systems. To account for additional pp interactions from the same and close-by bunch crossings, a set of minimum-bias interactions generated using Pythia  [28] 8.186 and the MSTW2008LO [29, 30] parton distribution function (PDF) set is superimposed onto the hard scattering events to reproduce the observed distribution of the average number of interactions per bunch crossing.

The signal samples are generated using MadGraph5_aMC@NLO  [31] v2.2.3 interfaced to Pythia  8.186 with the A14 [32] set of parameters (tune) for the modelling of the parton showering (PS), hadronization and underlying event. The matrix element (ME) calculation is performed at tree level and includes the emission of up to two additional partons. The PDF set used for the generation is NNPDF23LO [33]. The ME–PS matching is done using the CKKW-L [34] prescription, with a matching scale set to one quarter of the bottom squark mass. The cross-sections used to evaluate the signal yields are calculated to next-to-leading-order accuracy in the strong coupling constant, adding the resummation of soft gluon emission at next-to-leading-logarithmic accuracy (NLO + NLL) [3537].

SM background samples are simulated using different MC generator programs depending on the process. Events containing W or Z bosons with associated jets, including jets from the fragmentation of b- and c-quarks (heavy-flavour jets), are simulated using the Sherpa 2.1.1 [38] generator. Matrix elements are calculated for up to two additional partons at next-to-leading order (NLO) and four partons at leading order (LO) using the Comix  [39] and OpenLoops  [40] matrix element generators and merged with the Sherpa PS [41] using the ME + PS@NLO prescription [42]. The CT10 [43] PDF set is used in conjunction with a dedicated PS tune developed by the Sherpa authors. The W/Z+jets events are normalized to their next-to-NLO (NNLO) QCD theoretical cross-sections [44]. For a cross-check of the modelling of \(Z\)+jets background, samples of events containing a photon produced in association with jets, including heavy-flavour jets, are simulated using the Sherpa 2.1.1 generator with matrix elements calculated at LO with up to four partons for the LO calculation.

Diboson processes are also simulated using the Sherpa 2.1.1 generator with the same settings as the single-boson samples. They are calculated for up to one (ZZ) or zero (WWWZ) additional partons at NLO and up to three additional partons at LO. The NLO generator cross-sections are scaled down by 9 % to account for the usage of \(\alpha ^{}_\mathrm{QED}\)=1/129 rather than 1/132, the latter corresponding to the use of the parameters defined by the Particle Data Group as input to the \(G^{}_{\mu }\) scheme [45].

Top-quark pair and single-top-quark production (Wt- and s-channel) events are simulated using the Powheg-Box v2 [46] generator as described in Ref. [47] with the CT10 PDF set in the matrix element calculations. Electroweak t-channel single-top-quark events are generated using the Powheg-Box v1 generator. The top-quark mass is set to 172.5 GeV. For the \(t\bar{t}\) production, the \(h^{}_\mathrm{damp}\) parameter, which controls the transverse momentum (\(p_{\text {T}}\)) of the first additional emission beyond the Born configuration and thus regulates the \(p_{\text {T}}\)  of the recoil emission against the \(t\bar{t}\)  system, is set to the mass of the top quark. For all processes, the parton shower, fragmentation, and the underlying event are simulated using Pythia  6.428 [48] with the CTEQ6L1 PDF set and the corresponding Perugia 2012 tune (P2012) [49]. The \(t\bar{t}\) samples are normalized to their NNLO cross-section including the resummation of soft gluon emission at next-to-NLL accuracy using Top++2.0 [50]. Samples of single-top-quark production are normalized to the NLO cross-sections reported in Refs. [5153] for the s-, t- and Wt-channels, respectively. The associated production of a top-quark pair with a vector boson, \(t\bar{t} +W/Z \), is generated at LO with MadGraph5_aMC@NLO  v2.2.3 interfaced to Pythia  8.186, with up to two (\(t\bar{t} +W \)), one (\(t\bar{t} +Z \)) or no (\(t\bar{t} +W W \)) extra partons included in the matrix elements. The samples are normalized to their NLO cross-sections [31].

The EvtGen 1.2.0 program [54] is used for modelling the properties of bottom- and charm-hadron decays in all samples generated with MadGraph5_aMC@NLO and Powheg-Box.

Several samples produced without detector simulation are employed to derive systematic uncertainties associated with the specific configuration of the MC generators used for the nominal SM background samples. They include variations of the renormalization and factorization scales, the CKKW-L matching scale, as well as different PDF sets and fragmentation/hadronization models. Details of the MC modelling uncertainties are discussed in Sect. 7.

4 Event reconstruction

The search for bottom squark pair production is based on a selection of events with jets and large missing transverse momentum in the final state. Events containing electrons or muons are explicitly vetoed in the signal and validation regions, and are used to define control regions. An overlap removal procedure is applied to prevent double-counting of reconstructed objects. The details of the selections and overlap removal are given below.

Interaction vertices from the pp collisions are reconstructed from tracks in the inner detector. Events must have at least one primary reconstructed vertex, required to be consistent with the beamspot envelope and to consist of at least two tracks with \(p_{\text {T}}\) > 0.4 GeV. When more than one such vertex is found the one with the largest sum of the square of transverse momenta of associated tracks [55] is chosen.

Jet candidates are reconstructed from three-dimensional energy clusters [56] in the calorimeter using the anti-\(k_t\) jet algorithm [57] with a radius parameter of 0.4. Each topological cluster’s energy is calibrated to the electromagnetic scale prior to jet reconstruction. The reconstructed jets are then calibrated to the particle level by applying a jet energy scale (JES) derived from simulation and in situ corrections based on 8 TeV data [58, 59]. Quality criteria are imposed to identify jets arising from non-collision sources or detector noise and any event containing such a jet is removed [60]. Further track-based selections are applied to reject jets with \(p_{\text {T}} <60\) GeV and \(|\eta |<2.4\) that originate from pile-up interactions [61] and the expected event average energy contribution from pile-up clusters is subtracted using a factor dependent on the jet area [58]. Jets are classified as “baseline” and “signal”. Baseline jets are required to have \(p_{\text {T}} >20\) GeV and \(|\eta |<2.8\). Signal jets, selected after resolving overlaps with electrons and muons, are required to pass the stricter requirement of \(p_{\text {T}} >35\) GeV.

Jets are identified as \(b\)-jets if tagged by a multivariate algorithm which uses information about the impact parameters of inner detector tracks matched to the jet, the presence of displaced secondary vertices, and the reconstructed flight paths of b- and c-hadrons inside the jet [62]. The \(b\)-tagging working point with a 77 % efficiency, as determined in a simulated sample of \(t\bar{t}\) events, was chosen as part of the optimization process discussed in Sect. 5. The corresponding rejection factors against jets originating from \(c\)-quarks and from light quarks and gluons at this working point are 4.5 and 130, respectively [63]. To compensate for differences between data and MC simulation in the \(b\)-tagging efficiencies and mis-tag rates, correction factors derived from data-driven methods are applied to the simulated samples [64]. Candidate b-tagged jets are required to have \(p_{\text {T}} > 50\) GeV and \(|\eta |<\) 2.5.

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 [6567]. They are also required to lie within the fiducial volume \(|\eta |<2.47\). Muon candidates are reconstructed from matching tracks in the inner detector and the muon spectrometer. Events containing one or more muon candidates that have a transverse (longitudinal) impact parameter with respect to the primary vertex larger than 0.2 mm (1 mm) are rejected to suppress cosmic rays. Muon candidates are also required to satisfy “medium” quality criteria [68] and have \(|\eta |<\) 2.5. All electron and muon candidates must have \(p_{\text {T}}>\) 10 GeV. Lepton candidates remaining after resolving the overlap with baseline jets (see next paragraph) are called “baseline” leptons. In the control regions where lepton identification is required, “signal” leptons are chosen from the baseline set with \(p_{\text {T}}\) > 26 GeV to ensure full efficiency of the trigger and are required to be isolated from other activity using a selection designed to accept 99 % of leptons from Z boson decays. Signal electrons are further required to satisfy “tight” quality criteria [65]. Electrons (muons) 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 for both the electrons and muons. The MC events are corrected to account for differences in the lepton trigger, reconstruction and identification efficiencies between data and MC simulation.

The sequence to resolve overlapping objects begins by removing electron candidates sharing an inner detector track with a muon candidate. Next, jet candidates within \(\Delta R=\sqrt{(\Delta y)^2 + (\Delta \phi )^2} =0.2\) of an electron candidate are discarded, unless the jet is \(b\)-tagged, in which case the electron is discarded since it is likely to originate from a semileptonic \(b\)-hadron decay. Finally, any lepton candidate remaining within \(\Delta R = 0.4\) of any surviving jet candidate is discarded, except for the case where the lepton is a muon and the number of tracks associated with the jet is less than three.

The missing transverse momentum \({{\varvec{p}}}_{\mathrm {T}}^{\mathrm {miss}}\), with magnitude \(E_{\mathrm {T}}^{\mathrm {miss}}\), is defined as the negative vector sum of the \(p_{\text {T}}\) of all selected and calibrated physics objects 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. This soft term is calculated from inner detector tracks with \(p_{\text {T}}\) above 400 MeV matched to the primary vertex to make it more robust against pile-up contamination [69, 70].

Reconstructed photons, although not used in the main signal event selection, are selected in the regions employed in one of the alternative methods used to check the \(Z\)+jets background, as explained in Sect. 6. Photon candidates are required to have \(p_{\text {T}}\)  \( > 130\) GeV and \(|\eta | < 2.37\), to satisfy the tight photon shower shape and electron rejection criteria [71], and to be isolated.

5 Event selection

The selection of events is similar to that used in Ref. [72] and is based on the definition of one set of three overlapping signal regions (SRA) and a fourth distinct signal region (SRB). These were re-optimized for 3.2 fb\(^{-1}\) of 13 TeV pp collisions, and the selection criteria are summarized in Table 1. Only events with \(E_{\mathrm {T}}^{\mathrm {miss}}\) \(>250\) GeV are retained to ensure full efficiency of the trigger. Jets are ordered according to decreasing \(p_{\text {T}}\). Contamination from backgrounds with high jet multiplicity, particularly \(t\bar{t}\) production, is suppressed by vetoing events with a fourth jet with \(p_{\text {T}} >50\) GeV. To discriminate against multijet background, events where \(E_{\mathrm {T}}^{\mathrm {miss}}\) is aligned with a jet in the transverse plane are rejected by requiring \(\Delta \phi ^j_{\mathrm {min}}>0.4\), where \(\Delta \phi ^j_{\mathrm {min}}\) is the minimum azimuthal distance between \(E_{\mathrm {T}}^{\mathrm {miss}}\) and the leading four jets, and by requiring \(E_{\mathrm {T}}^{\mathrm {miss}}/m_{\mathrm {eff}}>0.25\), where \(m_{\mathrm {eff}}\) is defined as the scalar sum of the \(E_{\mathrm {T}}^{\mathrm {miss}} \) and the \(p_{\text {T}} \) of the two leading jets. For signal events, no isolated charged leptons are expected in the final state, and events where a baseline electron or muon is reconstructed are discarded.

Table 1 Summary of the event selection in each signal region. The term lepton is used in the table to refer to baseline electrons and muons. Jets (\(j_1\), \(j_2\), \(j_3\) and \(j_4\)) are labelled with an index corresponding to their decreasing order in \(p_{\text {T}}\)

The first set of signal regions, SRA, selects events with large \(E_{\mathrm {T}}^{\mathrm {miss}}\) where the two leading jets are \(b\)-tagged to target models with a large mass difference between the bottom squark and the neutralino. The main discriminating variable is the contransverse mass (\(m_\mathrm {CT}\)) [73], which is a kinematic variable that can be used to measure the masses of pair-produced semi-invisibly decaying heavy particles. For two identical decays of heavy particles into two visible particles (or particle aggregates) \(v_{1}\) and \(v_{2}\), and two invisible particles, \(m_\mathrm {CT} \) is defined as

$$\begin{aligned} m_\mathrm {CT} ^{2}(v_{1},v_{2}) = \left[ E_\mathrm{T}(v_{1}) + E_\mathrm{T}(v_{2}) \right] ^2 - \left[ {\mathbf{p }}_\mathrm{T}(v_{1}) - {\mathbf{p }}_\mathrm{T}(v_{2}) \right] ^2. \end{aligned}$$

In this analysis, \(v_1\) and \(v_2\) are the two leading \(b\)-jets. For signal events, these correspond to the \(b\)-jets from the squark decays and the invisible particles are the two neutralinos. The contransverse mass is invariant under equal and opposite boosts of the parent particles in the transverse plane. For systems of parent particles produced with small transverse boosts, \(m_\mathrm {CT}\) is bounded from above by an analytical combination of particle masses. This bound is saturated when the two visible objects are collinear. For \(t\bar{t}\) events, this kinematic bound is at 135 GeV, while for production of bottom squark pairs the bound is given by \(m_\mathrm{CT}^\mathrm{max} = (m^2_{\tilde{b}^{}_{1} } - m^2_{\tilde{\chi }^{0}_{1}} ) / m_{\tilde{b}^{}_{1} }\). The selection on \(m_\mathrm {CT}\) is optimized based on the bottom squark and neutralino masses considered and SRA is further divided into three overlapping regions, SRA250, SRA350 and SRA450, where the naming conventions reflects the minimum value allowed for \(m_\mathrm {CT}\). Finally, a selection on the invariant mass of the two \(b\)-jets (\(m_{bb}>200\) GeV) is applied to further enhance the signal yield over the SM background contributions. For a signal model corresponding to \(m^{}_{\tilde{b}^{}_{1} }\) = 800 GeV and \(m^{}_{\tilde{\chi }^{0}_{1}}\) = 1 GeV, 10, 8 and 5 % of the simulated signal events are retained by the SRA250, SRA350 and SRA450 selections, respectively.

The second type of signal region, SRB, selects events where a bottom squark pair is produced in association with a jet from initial-state radiation (ISR). The SRB targets models with a small mass difference between the \(\tilde{b}^{}_{1}\) and the \(\tilde{\chi }^{0}_{1}\), such that a boosted bottom squark pair is needed to satisfy the trigger requirements. Hence events are selected with large \(E_{\mathrm {T}}^{\mathrm {miss}}\), one high-\(p_{\text {T}}\) non-\(b\)-tagged leading jet and at least two additional \(b\)-jets. The leading jet is also required to be pointing in the direction opposite to the \(E_{\mathrm {T}}^{\mathrm {miss}}\) by requiring \(\Delta \phi (j_1,E_{\mathrm {T}}^{\mathrm {miss}}) > 2.5\), where \(\Delta \phi (j_1,E_{\mathrm {T}}^{\mathrm {miss}})\) is defined as the azimuthal angle between the leading jet and the \(E_{\mathrm {T}}^{\mathrm {miss}}\). For a signal model corresponding to \(m^{}_{\tilde{b}^{}_{1} }\) = 400 GeV and \(m^{}_{\tilde{\chi }^{0}_{1}}\) = 300 GeV, about 0.3 % of the simulated events are retained by the SRB selection.

Table 2 Definition of the control regions associated with SRA and SRB. Control regions are defined to study the contribution from \(Z\)+hf, \(t\bar{t}\), single top-quark and \(W\)+hf production in SRA and the contribution from \(Z\)+hf and \(t\bar{t}\) in SRB. The term lepton is used in the table to refer to signal electrons and muons. Jets (\(j_1\), \(j_2\), \(j_3\) and \(j_4\)) are labelled with an index corresponding to their decreasing order in \(p_{\text {T}}\). SFOS indicates the same-flavour opposite-sign two-lepton selection

6 Background estimation

The dominant SM background processes in the signal regions are the production of W or Z bosons in association with heavy-flavour jets (referred to as \(W\)+hf and \(Z\)+hf) and the production of top quarks. In particular, events with \(Z\)+hf production followed by the \(Z\rightarrow \nu \bar{\nu }\) decay have the same signature as the signal and are the largest (irreducible) background in SRA. The background in SRB is dominated by top-quark production in events with a charged lepton in the final state that is not reconstructed, either because the lepton is a hadronically decaying \(\tau \), or because the electron or muon is not identified or out of detector acceptance.

Monte Carlo simulation is used to estimate the background yield in the signal regions, after normalizing the Monte Carlo prediction for the major backgrounds to data in control regions (CR) constructed to enhance a particular background and to be kinematically similar but orthogonal to the signal regions. The control regions are defined by explicitly requiring the presence of one or two leptons (electrons or muons) in the final state and applying further selection criteria similar to those of the corresponding signal regions. To select events with good-quality electrons and muons the “signal lepton” selection is applied to them. Events with additional baseline lepton candidates are vetoed. For SRA, the normalizations of the top-quark pair, single-top-quark, \(Z\)+jets and \(W\)+jets backgrounds are estimated simultaneously by making use of four control regions, while for SRB two control regions are used to determine the normalization of the top-quark pair and \(Z\)+jets backgrounds since other contributions are sub-dominant. A likelihood function is built as the product of Poisson probability functions, using as constraints the observed and expected (from MC simulation) event yields in the control regions but not the yields in the corresponding SR [74]. The normalization factors for each of the simulated backgrounds are adjusted simultaneously via a profile likelihood fit [75] (referred to as the “background-only fit”). All systematic uncertainties (discussed in Sect. 7) are treated as nuisance parameters in the fit. The background normalization parameters are applied to the signal regions also taking into account correlations in the yield predictions between different regions.

Two same-flavour opposite-sign (SFOS) two-lepton (electron or muon) control regions with dilepton invariant mass near the Z boson mass (\(76< m_{\ell \ell } <106\) GeV) and two \(b\)-tagged jets provide data samples dominated by Z boson production. For these control regions, labelled in the following as CRzA and CRzB for SRA and SRB respectively, the \(p_{\text {T}}\) of the leptons is added vectorially to the \({{\varvec{p}}}_{\mathrm {T}}^{\mathrm {miss}}\) to mimic the expected missing transverse momentum spectrum of \(Z\rightarrow \nu \bar{\nu }\) events, and is indicated in the following as \(E_{\mathrm {T}}^{\mathrm {miss,cor}} \) (lepton corrected). In addition, the uncorrected \(E_{\mathrm {T}}^{\mathrm {miss}}\) of the event is required to be less than 100 (70) GeV in CRzA (CRzB) in order to further enhance the Z boson contribution. In the case of CRzA, a \(m_{bb}>200\) GeV selection is also imposed. Two control regions, CRttA and CRttB defined for SRA and SRB respectively, dominated by \(t\bar{t}\) production are identified by selecting events with exactly one lepton (\(e,\mu \)) and a set of requirements similar to those for SRA and SRB, with the additional requirement \(m_{bb}<200\) GeV in CRttA to separate the \(t\bar{t}\) pair contribution from single top-quark production. To assist in the estimation of the background from \(W\)+hf and single top-quark production in SRA two further one-lepton control regions (CRwA and CRstA) are defined. These control regions exploit kinematic features to differentiate these processes from \(t\bar{t}\) production. For CRstA a selection is applied to the minimum invariant mass of either of the \(b\)-jets and the charged lepton (\(m_{b\ell }^{\mathrm {min}}\)). For CRwA only one \(b\)-jet is required, the selection on \(m_{bb}\) is replaced by a selection on the invariant mass of the two leading jets (\(m_{bj}\)) and a further selection is applied to the event transverse mass \(m_{\mathrm {T}}\), defined as \(m_{\mathrm {T}} = \sqrt{2 {p}_\mathrm {T}^\mathrm{lep} E_{\mathrm {T}}^{\mathrm {miss}}- 2 {{\varvec{p}}}_{\mathrm {T}}^{\mathrm {lep}} \cdot {{\varvec{p}}}_{\mathrm {T}}^{\mathrm {miss}}} \). The definitions of the control regions are summarized in Table 2.

The contributions from diboson and \(t\bar{t} + W/Z\) processes are minor and they are collectively called “Others” in the following. They are estimated from MC simulation for both the signal and the control regions and included in the fit procedure, and are allowed to vary within their uncertainty. The background from multijet production is estimated from data using a procedure described in detail in Ref. [76] and modified to account for the flavour of the jets. The procedure consists of smearing the jet response in events with well-measured \(E_{\mathrm {T}}^{\mathrm {miss}}\) (seed events). The jet response function is obtained from MC dijet events and cross-checked in events from data where the \(E_{\mathrm {T}}^{\mathrm {miss}}\) can be unambiguously attributed to the mis-measurement of one of the jets. The contribution from multijet production in all regions is found to be negligible.

Table 3 Fit results in the control regions associated with the SRA and SRB selection for an integrated luminosity of 3.2 fb\(^{-1}\). The results are obtained from the control regions using the background-only fit. The uncertainties include statistical, detector-related and theoretical systematic components. The individual uncertainties can be correlated and do not necessarily add in quadrature to the total systematic uncertainty. The pure MC estimate is used for backgrounds for which a dedicated CR is not defined, e.g. for smaller backgrounds (indicated as “Other”) and for single top-quark and \(W\)+jets production in CRzB and CRttB. A dash indicates a negligible background

The results of the background-only fit are shown in Table 3, where the contribution from individual backgrounds is shown separately as a purely MC-based prediction and with the rescaling from the fit procedure. The \(m_{bb}\) distribution in CRzA and the \(p_{\text {T}}\) of the leading jet in CRttB are shown in Fig. 1 after the backgrounds were rescaled as a result of the background-only fit, showing good agreement in the shape of the distributions in the two control regions used to estimate the dominant backgrounds in SRA and SRB.

Fig. 1
figure 1

Left \(m_{bb}\) distribution in CRzA before the final selection \(m_{bb}>200\) GeV is applied (indicated by the arrow). Right \(p_{\text {T}}\) of the leading jet in CRttB. In both distributions the MC normalization is rescaled using the results from the background-only fit, showing good agreement between data and the predicted SM shapes. The shaded band includes statistical and detector-related systematic uncertainties as detailed in Sect. 7 and the last bin includes overflows

The full background estimation procedure is validated by comparing the background predictions and the shapes of the distributions of the key analysis variables from the fit results to those observed in dedicated validation regions. They are defined to be mutually exclusive and kinematically similar to the signal regions, with low potential contamination from signal. For SRA, two validation regions are defined by using the same definition as SRA250, and inverting the selections on \(m_{bb}\) and \(m_\mathrm {CT}\). The SM predictions are found to overestimate the data by about one standard deviation. For SRB, a validation region is defined by selecting events with 250 \(<E_{\mathrm {T}}^{\mathrm {miss}}<\) 300 GeV and good agreement is observed between data and predictions.

As a further validation, two alternative methods are used to estimate the \(Z\)+hf contribution. The first method exploits the similarity of the \(Z\)+jets and \(\gamma +\)jets processes [76]. For \(p_{\text {T}}\) of the photon significantly larger than the mass of the Z boson, the kinematics of \(\gamma +\)jets events strongly resemble those of \(Z\)+jets events. The event yields are measured in control regions identical to the SRA and SRB, with the \(E_{\mathrm {T}}^{\mathrm {miss}}\)-based selections replaced by selections on the \(p_{\text {T}}\) of the photon vectorially added to the \({{\varvec{p}}}_{\mathrm {T}}^{\mathrm {miss}}\). The yields are then propagated to the actual SRA and SRB using a reweighting factor derived using the MC simulation. This factor takes into account the different kinematics of the two processes and residual effects arising from the acceptance and reconstruction efficiency for photons.

In the second alternative method, applied to SRA only where the \(Z\)+hf contribution is dominant, the MC simulation is used to verify that the shape of the \(m_\mathrm {CT}\) distribution for events with no \(b\)-tagged jets is compatible with the shape of the \(m_\mathrm {CT}\) distribution for events where two b-tagged jets are present. A new highly populated \(Z\)+jets CR is defined selecting \(Z\rightarrow \ell \ell \) events with no b-tagged jets. The \(m_\mathrm {CT}\) distribution in this CR is constructed using the two leading jets and is used to estimate the shape of the \(m_\mathrm {CT}\) distribution in the SRA, whilst the normalization in SRA is rescaled based on the ratio in data of \(Z \rightarrow \ell \ell \) events with no \(b\)-tagged jets to events with two \(b\)-tagged jets. Additional MC-based corrections are applied to take into account the two-lepton selection in this CR.

Table 4 Estimated \(Z\)+jets yields in the signal regions as obtained using the default and the two alternative methods. The errors include all the uncertainty sources discussed in Sect. 7. The “MC-based post-fit” uncertainty does not include the additional uncertainty to account for the difference between the three methods

The two alternative methods are in agreement within uncertainties with the estimates obtained with the profile likelihood fit to the control regions (Table 4). Experimental and theoretical systematic uncertainties in the nominal and alternative method estimates are taken into account (see Sect. 7). The difference between the alternative methods and the background-only fit is taken into account as an additional systematic uncertainty in the final \(Z\)+hf yields.

7 Systematic uncertainties

Several sources of experimental and theoretical systematic uncertainty are considered in this analysis. Their impact is reduced through the normalization of the dominant backgrounds in the control regions with kinematic selections resembling those of the corresponding signal region (see Sect. 6). Uncertainties due to the limited number of events in the CRs are also taken into account in the fit. The individual contributions are outlined in Table 5.

Table 5 Summary of the dominant experimental and theoretical uncertainties for each signal region. Uncertainties are quoted as relative to the total uncertainty, with a range indicated for the three SRAs. For theoretical modelling, uncertainties per dominant SM background process are quoted. The individual uncertainties can be correlated, and do not necessarily add in quadrature to the total background uncertainty

The dominant detector-related systematic effects are due to the uncertainties in the jet energy scale (JES) [58] and resolution (JER) [59], and in the b-tagging efficiency and mis-tagging rates. The JES and JER uncertainties are estimated from 13 TeV data, while the uncertainties related to \(b\)-tagging are estimated from 8 TeV data and extrapolated to 13 TeV and to the Run-2 detector conditions. The uncertainties associated with lepton and photon reconstruction and energy measurements are also considered but have a small impact on the final results. Lepton, photon and jet-related uncertainties are propagated to the calculation of the \(E_{\mathrm {T}}^{\mathrm {miss}}\), and additional uncertainties are included in the energy scale and resolution of the soft term. The overall experimental uncertainty in the SM background estimate is found to be around 20 % for the SRAs and 15 % for the SRB.

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 uncertainty arises from Z+jets MC modelling for SRA and \(t\bar{t} \) modelling for SRB. The Z+jets (as well as W+jets) modelling uncertainties are evaluated using alternative samples generated with different renormalization and factorization scales, merging (CKKW-L) and resummation scales. An additional one-sided uncertainty in the \(Z\)+jets estimate is taken as the largest deviation between the nominal background-only fit result and each of the alternative data-driven estimates described in Sect. 6. This results in an additional 25, 25 and 40 % one-sided uncertainty in SRA250, SRA350 and SRA450, respectively. Finally, a 40 % uncertainty [77] is assigned to the heavy-flavour jet content in W+jets, estimated from MC simulation in SRB. For SRA, the uncertainty accounts for the different requirements on \(b\)-jets between CRwA and the signal region.

Table 6 Fit results in all signal regions for an integrated luminosity of 3.2 fb\(^{-1}\). The results are obtained from the control regions. The background normalization parameters obtained with the background-only fit are applied to the SRs. The individual uncertainties, including detector-related and theoretical systematic components, are symmetrized and can be correlated and do not necessarily add in quadrature to the total systematic uncertainty

Uncertainties in the modelling of the top-quark pair and single-top-quark (Wt) backgrounds are sub-dominant in SRA and dominant in SRB. They are computed as the difference between the predictions from nominal samples and those of additional samples differing in generator or parameter settings. Hadronization and PS uncertainties are estimated using samples generated with Powheg-Box v2 and showered by Herwig++ v2.7.1 [78] with the UEEE5 underlying event tune. Uncertainties related to initial- and final-state radiation modelling, tune and (for \(t\bar{t}\)  only) choice of \(h^{}_\mathrm{damp}\) parameter in Powheg-Box v2 are evaluated using alternative settings of the generators. Finally, an alternative generator MadGraph5_aMC@NLO with showering by Herwig++ v2.7.1 is used to estimate the generator uncertainties. Uncertainties in smaller backgrounds such as diboson and ttV are also estimated by comparisons of the nominal sample with alternative samples differing in generator or parameter settings (Powheg v2 with showering by Pythia  8.210 for diboson, renormalization and factorization scale and A14 tune variations for ttV) and are found to be negligible. The cross-sections used to normalize the MC yields to the highest order available are varied according to the scale uncertainty of the theoretical calculation, i.e. 5 % for W, Z boson and top-quark pair production, 6 % for diboson, 13 and 12 % for ttW and ttZ, respectively.

For the SUSY signal processes, both the experimental and theoretical uncertainties in the expected signal yield are considered. Experimental uncertainties, which are found to be between 20 and 25 % across the \(\tilde{b}^{}_{1} \)\(\tilde{\chi }^{0}_{1}\) mass plane for all SRs, are largely dominated by the uncertainty in the b-tagging efficiency in SRA, while JER and b-tagging uncertainties are dominant in SRB with equal contributions. Theoretical uncertainties in the NLO + NLL cross-section are calculated for each SUSY signal scenario and are dominated by the uncertainties in the renormalization and factorization scales, followed by the uncertainty in the PDF. They vary between 15 and 20 % for bottom squark masses in the range between 400 and 900 GeV. Additional uncertainties in the modelling of initial-state radiation in SUSY signal MC samples are taken into account and contribute up to 5 %.

8 Results and interpretation

Table 6 reports the observed number of events and the SM predictions after the background-only fit for each signal region. The largest background contribution in SRA arises from \(Z\rightarrow \nu \bar{\nu }\) produced in association with b-quarks whilst top-quark pair production dominates SM predictions for SRB. The background-only fit results are compared to the pre-fit predictions based on MC simulation. Figures 2 and 3 show the comparison between the observed data and the SM predictions for some relevant kinematic distributions in SRA250 and SRB, respectively, prior to the selection on the variable shown. For illustrative purposes, the distributions expected for a scenario with bottom squark and neutralino masses of 700 GeV (400 GeV) and 1 GeV (300 GeV), respectively, are shown for SRA250 (SRB). No excess above the expected Standard Model background yield is observed. The results are translated into upper limits on contributions from new physics beyond the SM (BSM) for each signal region. The CL\(_s\) method [79, 80] is used to derive the confidence level of the exclusion; signal models with a CL\(_s\) value below 0.05 are said to be excluded at 95 % CL. The profile-likelihood-ratio test statistic is used to exclude the signal-plus-background hypothesis for specific signal models. When normalized by the integrated luminosity of the data sample, results can be interpreted as corresponding upper limits on the visible cross-section, \(\sigma ^{}_\mathrm{vis}\), defined as the product of the BSM production cross-section, the acceptance and the selection efficiency of a BSM signal. Table 7 summarizes 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 ^{}_\mathrm{vis}\).

Fig. 2
figure 2

Left \(m_\mathrm {CT}\) distribution in SRA250 with all the selection criteria applied except the \(m_\mathrm {CT}\) threshold. Right \(m_{bb}\)  distribution in SRA250 with all selection criteria applied except the \(m_{bb}\) requirement. The arrows indicate the final applied selection. The shaded band includes statistical and detector-related systematic uncertainties. The SM backgrounds are normalized to the values determined in the fit. For illustration the distributions expected for one signal model with bottom squark and neutralino masses of 700 and 1 GeV, respectively, are overlaid. The last bin includes overflows

Fig. 3
figure 3

Left \(E_{\mathrm {T}}^{\mathrm {miss}}\) distribution in SRB with all the selection criteria applied except the \(E_{\mathrm {T}}^{\mathrm {miss}}\) threshold. Right leading jet \(p_{\text {T}}\) distribution in SRB with all the selection criteria applied except the selection on the leading jet \(p_{\text {T}}\) itself. The arrows indicate the final selection applied in SRB. The shaded band includes statistical and detector-related systematic uncertainties. The SM backgrounds are normalized to the values determined in the fit. For illustration the distributions expected for one signal model with bottom squark and neutralino masses of 400 and 300 GeV, respectively, are overlaid. The last bin includes overflows

Exclusion limits are obtained assuming a specific SUSY particle mass hierarchy such that the lightest bottom squark decays exclusively via \(\tilde{b}^{}_{1} \rightarrow b \tilde{\chi }^{0}_{1}\). In this case, the fit procedure takes into account not only correlations in the yield predictions between control and signal regions due to common background normalization parameters and systematic uncertainties but also contamination in the control regions from SUSY signal events (found to be generally negligible). The experimental systematic uncertainties in the signal are taken into account for this calculation and are assumed to be fully correlated with those in the SM background. Figure 4 shows the observed (solid line) and expected (dashed line) exclusion contours at 95 % CL in the \(\tilde{b}^{}_{1} \)\(\tilde{\chi }^{0}_{1}\) mass plane.

Table 7 Left to right: 95 % CL upper limits on the visible cross-section (\(\langle \epsilon A \sigma \rangle _\mathrm{obs}^{95}\)) and on the number of signal events (\(S_\mathrm{obs}^{95}\) ). The third column (\(S_\mathrm{exp}^{95}\)) shows the 95 % CL upper limit on the number of signal events, given the expected number (and \(\pm 1\sigma \) variations of the expectation) of background events
Fig. 4
figure 4

Observed and expected exclusion limits at 95 % CL, as well as \(\pm 1 \sigma \) variation of the expected limit, in the \(\tilde{b}^{}_{1} \)\(\tilde{\chi }^{0}_{1}\) mass plane. The SR with the best expected sensitivity is adopted for each point of the parameter space. The yellow band around the expected limit (dashed line) shows the impact of the experimental and SM background theoretical uncertainties. The dotted lines show the impact on the observed limit of the variation of the nominal signal cross-section by \(\pm 1 \sigma \) of its theoretical uncertainties. The exclusion limits from the Run-1 ATLAS searches [72, 81] and from the 13 TeV monojet search [82] are also superimposed. The latter limit is only published for values of \(m_{\tilde{b}^{}_{1} }- m_{\tilde{\chi }^{0}_{1}}=\) 5 and 20 GeV

At each point of the parameter space, the SR with the best expected sensitivity is adopted. Sensitivity to scenarios with the largest mass difference between the \(\tilde{b}^{}_{1} \) and the \(\tilde{\chi }^{0}_{1}\) is achieved with the most stringent \(m_\mathrm {CT}\)  threshold (SRA450). Sensitivity to scenarios with smaller mass differences is achieved predominantly with searches based on the presence of a high-\(p_{\text {T}}\) ISR jet, as in the dedicated SRB selection and the search described in Ref. [81]. Bottom squark masses up to 800 (840) GeV are excluded for \(\tilde{\chi }^{0}_{1}\) masses below 360 (100) GeV. Differences in mass above 100 GeV between \(\tilde{b}^{}_{1} \) and \(\tilde{\chi }^{0}_{1}\) are excluded up to \(\tilde{b}^{}_{1} \) masses of 500 GeV. The expected exclusion constraints are about 20 GeV lower than the observation for high bottom squark masses and about 60 GeV lower for SUSY models with small \(\tilde{b}^{}_{1} \)\(\tilde{\chi }^{0}_{1}\) mass splitting. The current results significantly extend the \(\sqrt{s} = 8\) TeV limits [72, 81] despite the lower integrated luminosity mostly because of the increase in centre-of-mass energy of the LHC. Furthermore, the sensitivity of the analysis benefits from the advanced algorithms adopted to identify b-jets and use of information from the newly installed IBL pixel layer in the Run-2 ATLAS detector, as well as from improved techniques to estimate SM background contributions and their systematic uncertainties.

9 Conclusion

In summary, the results of a search for bottom squark pair production are reported. The analysis uses 3.2 fb\(^{-1}\)  of pp collisions at \(\sqrt{s}=13\) TeV collected by the ATLAS experiment at the Large Hadron Collider in 2015. Bottom squarks are searched for in events containing large missing transverse momentum and up to four jets, exactly two of which are identified as \(b\)-jets. No excess above the expected Standard Model background yield is found. Exclusion limits at 95 % confidence level are placed on the visible cross-section and on the mass of the bottom squark in phenomenological supersymmetric R-parity-conserving models in which the \(\tilde{b}^{}_{1} \) is the lightest squark and is assumed to decay exclusively via \(\tilde{b}^{}_{1} \rightarrow b \tilde{\chi }^{0}_{1}\), where \(\tilde{\chi }^{0}_{1}\) is the lightest neutralino. Bottom squark masses up to 800 GeV are excluded for \(\tilde{\chi }^{0}_{1}\) masses below 360 GeV (840 GeV for \(\tilde{\chi }^{0}_{1}\) masses below 100 GeV) whilst differences in mass above 100 GeV between \(\tilde{b}^{}_{1} \) and \(\tilde{\chi }^{0}_{1}\) are excluded up to \(\tilde{b}^{}_{1} \) masses of 500 GeV. The results significantly extend the constraints on bottom squark masses with respect to Run-1 searches.