Search for dark matter produced in association with a Standard Model Higgs boson decaying into b-quarks using the full Run 2 dataset from the ATLAS detector

The production of dark matter in association with Higgs bosons is predicted in several extensions of the Standard Model. An exploration of such scenarios is presented, considering final states with missing transverse momentum and b-tagged jets consistent with a Higgs boson. The analysis uses proton-proton collision data at a centre-of-mass energy of 13 TeV recorded by the ATLAS experiment at the LHC during Run 2, amounting to an integrated luminosity of 139 fb−1. The analysis, when compared with previous searches, benefits from a larger dataset, but also has further improvements providing sensitivity to a wider spectrum of signal scenarios. These improvements include both an optimised event selection and advances in the object identification, such as the use of the likelihood-based significance of the missing transverse momentum and variable-radius track-jets. No significant deviation from Standard Model expectations is observed. Limits are set, at 95% confidence level, in two benchmark models with two Higgs doublets extended by either a heavy vector boson Z′ or a pseudoscalar singlet a and which both provide a dark matter candidate χ. In the case of the two-Higgs-doublet model with an additional vector boson Z′, the observed limits extend up to a Z′ mass of 3 TeV for a mass of 100 GeV for the dark matter candidate. The two-Higgs-doublet model with a dark matter particle mass of 10 GeV and an additional pseudoscalar a is excluded for masses of the a up to 520 GeV and 240 GeV for tan β = 1 and tan β = 10 respectively. Limits on the visible cross-sections are set and range from to 0.05 fb to 3.26 fb, depending on the missing transverse momentum and b-quark jet multiplicity requirements.


Introduction
Various astrophysical observations based on gravitational interactions [1,2] strongly support the existence of dark matter (DM) which interacts through neither the strong nor the electromagnetic force. However, the Standard Model of particle physics (SM) provides no suitable DM candidate particle. There are many complementary search strategies for DM including direct-detection [3][4][5][6][7] and indirect-detection experiments [8] as well as searches at particle colliders.
Since DM particles do not interact electromagnetically or strongly even direct-detection experiments have low efficiencies. Therefore, instead of attempting to detect them directly, any DM particles produced in proton-proton collisions at the Large Hadron Collider [9] (LHC) would be deduced from an imbalance in the transverse momentum measured in that collision event (E miss T ). This means that DM particles can only be detected if they are produced in association with visible particles. When there is only one such particle, this gives rise to event topologies referred to as 'mono-X' final states, where X refers Mono-X topologies are typically dominated by cases where the visible particle is produced from initial-state radiation (ISR), as the couplings between the SM particle X and the initial-state quarks or gluons are much larger than its couplings to the final-state DM particles. A counterexample is the case where the visible particle is a SM Higgs boson (named the mono-Higgs signature). Given that the coupling of the Higgs boson to light quarks and gluons is highly suppressed, a Higgs boson is more likely to be produced through final-state radiation (FSR) or as part of the same process that produces the DM particles. This means that this topology is only sensitive to models where the Higgs boson couples directly to DM or some other beyond-the-SM (BSM) particle involved in DM production. However, in these cases the DM-SM interaction is probed directly [19], potentially providing more information about the structure of the DM-SM coupling in the event of a discovery. There are also many models where the coupling between BSM particles and the Higgs boson is enhanced, for instance models where DM is connected to electroweak symmetry breaking [20,21] or where DM particles couple to the SM only through the Higgs sector (Higgs portal models) [22]. These features make the mono-Higgs signature an important part of the LHC DM search programme.
The two-Higgs-doublet model (2HDM) [23] extends the SM with a second Higgs doublet. This predicts a total of five Higgs bosons after mixing: two charged scalars H ± , two neutral CP-even scalars h and H, and one neutral CP-odd scalar A. Two simplified benchmark signal models are used: the Z -2HDM [24] and 2HDM+a [25,26]. In all models considered here the mass of the lighter neutral CP-even scalar h is required to match that of the Higgs boson observed at the LHC and the Yukawa couplings are defined according to the Type II 2HDM.
The Z -2HDM has an additional heavy vector boson, denoted Z , whose coupling to quarks g Z and mass m Z are free parameters. The production mechanism for the mono-Higgs signature is shown in figure 1. DM is introduced into the model as a new fermion which couples to the CP-odd scalar A with a coupling strength denoted by g χ . The coupling strength and the DM particle mass m χ are treated as free parameters. This model is used mainly as a benchmark for high-mass resonances.
The 2HDM+a scenario is the simplest renormalisable and gauge-invariant extension of a simplified pseudoscalar mediator model. It adds a new pseudoscalar singlet which mediates the interactions between the SM and a singlet fermion χ identified as the DM candidate. The coupling between the pseudoscalar and DM singlets, denoted y χ , and the  mass of the DM m χ are free parameters. This singlet mixes with the pseudoscalar A from the two Higgs doublets, with the mixing angle θ and the mass of the resulting pseudoscalar a being free parameters of the model. A major advantage of the 2HDM+a scenario over simpler models is that it generates a wider variety of experimental signatures which can provide complementary exclusion regions from different types of experiment. There are two main production mechanisms for the mono-Higgs signature in this model, as shown in figure 2.
In the type-II 2HDM considered in this paper the coupling between down-type quarks and the A boson scales with tan β, the ratio of the vacuum expectation values of the two Higgs doublets. This means that for low tan β values (tan β 5) the gluon-gluon fusion (ggF) mechanism shown in figure 2(a) dominates, whereas for higher tan β values the bassociated production (bbA) shown in figure 2(b) is dominant. Signal grids are generated where each of the two production mechanisms are used exclusively. For each grid a tan β value is chosen that ensures that the corresponding production mechanism is dominant: tan β = 1 for the ggF grid and tan β = 10 for the bbA grid.
Similar analyses were performed using data taken during the years 2015-2016 by AT-LAS [27] and CMS [28]. Other mono-Higgs analyses were also performed on the same datasets in final states where the Higgs boson decays into a pair of photons in ATLAS [29] or either a pair of photons, a τ + τ − pair in CMS [30] or a pair of W or Z bosons [28]. Beyond the large increase in integrated luminosity (from 36 fb −1 to 139 fb −1 ) a number of analysis improvements extend the sensitivity beyond the previous ATLAS search.
In the previous ATLAS search, events were required to have either one or two b-jets, whereas in this search the events are divided into regions with either exactly two or at least three b-jets. Introducing the exclusive three b-jet category improves the sensitivity to the 2HDM+a bbA production mechanism which was not considered in the previous ATLAS search, while the one b-jet category does not provide any significant improvement due to large backgrounds with high uncertainties.
The analysis also benefits from using particle-flow objects for jet reconstruction [31], neural-network based b-jet [32] and τ -lepton [33] identification, and variable-radius trackjets [34] to identify boosted Higgs boson candidates. In order to reduce backgrounds containing fake E miss JHEP11(2021)209

Data and simulated event samples
This search uses 139 fb −1 of proton-proton collision data recorded by the ATLAS detector at a centre-of-mass energy of 13 TeV during the years 2015-2018 (Run 2). The uncertainty in the total integrated luminosity for the full Run 2 dataset is 1.7% [41], obtained using the LUCID-2 detector [42] for the primary luminosity measurements. All events used are required to pass basic data-quality requirements which ensure that all components of the ATLAS detector were functioning correctly [43]. Events selected for the analysis search regions were collected by the primary E miss T triggers [44], which select those that have a large transverse momentum imbalance. The algorithms used to calculate the imbalance, and the thresholds required, varied during the data-taking. The 'primary' E miss T trigger in a run is the most efficient available E miss T trigger, in all cases reaching full efficiency by an offline E miss T value of approximately 200 GeV. Simulated event samples corresponding to the Z -2HDM signal were generated at leading order (LO) in QCD in the 5-flavour scheme using Mad-Graph5_aMC@NLO v2.6.5 [45] interfaced to Pythia 8.240 [46] using a set of tuned parameters called the A14 tune [47]. Following the recommendations of the ATLAS-CMS Dark Matter Forum [48] the coupling of the Z boson to quarks was fixed to g Z = 0.8, the mass of the DM candidate was set to m χ = 100 GeV, the Aχχ coupling g χ was set to 1, tan β was set to 1, and the alignment limit, i.e. sin(β − α) = 1, was assumed, where α is the mixing angle between the two CP-even Higgs bosons. The samples were generated with varying m Z between 600 and 3600 GeV and m A between 300 and 1300 GeV.
Simulated event samples corresponding to the 2HDM+a signal were generated at LO using MadGraph5_aMC@NLO v2.6.7 interfaced to Pythia 8.244 with the A14 tune. Samples were generated separately for the ggF and bbA production modes, with the former generated in the 4-flavour scheme setting tan β = 1 and the latter in the 5-flavour scheme setting tan β = 10. Following the recommendations of the LHC Dark Matter Working Group [25], the mass of the DM candidate was set to m χ = 10 GeV, the Yukawa coupling between the DM candidate and the pseudoscalar a was set to y χ = 1 and the Higgs quartic couplings were set to λ 3 = λ P 1 = λ P 2 = 3. The chosen value of m χ ensures that the a → χχ branching ratio is significant for all values of m a used. The pseudoscalar mixing angle was set to sin θ = 0.35 and the alignment limit was assumed. The samples were generated with varying m A between 250 and 2000 GeV and m a between 100 and 600 GeV.
For all signal samples the masses of the heavy Higgs bosons were considered degenerate (m A = m H = m H ± ). The mass of the lightest CP-even Higgs boson was set to match that of the Higgs boson discovered at the LHC, i.e. m h = 125 GeV [49].
Background events from the production of a single weak vector boson (V = W, Z) in association with jets or of a pair of weak bosons (diboson) were simulated using Sherpa 2.2.1 [50], with Sherpa 2.2.2 used for gg-initiated diboson production. For the V +jets samples, next-to-leading-order (NLO) matrix elements for up to two jets and LO matrix elements for up to four jets were calculated using the Comix [51] and Open-Loops [52,53] libraries. For the qq-initiated diboson samples, NLO matrix elements for up to one additional jet and LO matrix elements for up to three additional jets were used, -5 -

JHEP11(2021)209
while the gg-initiated processes were generated using LO matrix elements for up to one additional jet. The samples were matched with the Sherpa parton shower [54] using the MEPS@NLO prescription [55][56][57][58].
The W /Z+h samples were generated using PowhegBox v2. The cross-sections of the qq-initiated processes were calculated at NNLO QCD and NLO EW accuracy, using the Powheg MiNLO procedure [78,79]. The cross-sections of the gg-initiated processes were calculated at NLO+next-to-leading-logarithm (NLL) accuracy in QCD [80][81][82].
The ttV samples were generated using MadGraph5_aMC@NLO v2.3.3 at NLO. The cross-sections were calculated at NLO QCD and EW accuracies as provided by ref. [76].
All samples were generated using the NNPDF3.0NLO parton distribution function (PDF) set [83] apart from the Sherpa W /Z+jets and diboson samples, which were generated using the NNPDF3.0NNLO PDF set along with a dedicated tune developed by the Sherpa authors, and the t-channel single-top-quark production samples, which were generated using the NNPDF3.0NLONF4 PDF set. The PowhegBox samples were interfaced with Pythia 8.230 for the parton shower and hadronisation. Of these, the tt, single-top-quark and tth samples used the NNPDF2.3LO PDF set [83] and the A14 tune [84], while the W /Z+h samples used the CTEQ6L1 PDF set [85] and the AZNLO tune [86]. For the PowhegBox top-quark samples the decays of b-and c-hadrons were simulated using EvtGen 1.6.0 [87]. The ttV samples were interfaced with Pythia 8.210 using the same tune and PDF set as the PowhegBox tt, single-top-quark and tth samples, and using EvtGen 1.2.0 for the decays of b-and c-hadrons.
In order to simulate the effect of additional pp collisions in the same and neighbouring bunch crossings (pile-up) all samples were overlaid with multiple pp collisions simulated with Pythia 8.186 using the NNPDF2.3LO PDF set and the A3 tune [88]. The response of the detector was modelled with a detector simulation [89] based on Geant4 [90].

Object definitions
Primary vertex. Primary vertices are constructed using at least two ID tracks with p T > 500 MeV [91]. The primary vertex with the largest sum of squared track transverse momenta ( p 2 T ) is selected as the hard-scatter vertex, henceforth only referred to as the primary vertex.
Jets. Jets are reconstructed using the anti-k t algorithm [92,93]. The analysis considers three types of jets to better match the different event topologies. Small-radius (small--6 -JHEP11(2021)209 R) jets are constructed from particle-flow objects formed from ID tracks and calorimeter energy clusters [31] using a radius parameter of R = 0.4. This radius parameter is designed to capture jets initiated by a gluon, light quark or b-quark. Small-R jets are classified as central (|η| < 2.5) or forward (2.5 < |η| < 4.5). Central small-R jets are required to have p T > 20 GeV and forward small-R jets p T > 30 GeV. In order to remove the impact of jets predominantly formed from particles from pile-up vertices, central small-R jets are required to pass the 'Tight' jet vertex tagger (JVT) [94] working point (WP). 2 In topologies where the Higgs boson decay h → bb cannot be resolved into two small-R jets, large-radius (large-R) jets with a radius parameter of R = 1.0 are used, constructed from calorimeter energy clusters calibrated using the local hadronic cell weighting (LCW) scheme [95]. This radius parameter is chosen so that a single large-R jet should capture all jets produced in the decay of a boosted heavy object, such as a Higgs boson. To reduce the impact of pile-up, these jets are then 'trimmed', removing any R = 0.2 subjets which have less than 5% of the original jet energy [96]. In order to identify subjets originating from b-hadrons within the large-R jets, jets are also constructed from ID tracks, using a variant of the anti-k t algorithm with a radius parameter that shrinks as the p T of the proto-jet increases [34]. These are referred to as variable-radius (variable-R) track-jets and are matched to the large-R jets by ghost association [97]. The radius parameter is set to R = 30 GeV/p T , with minimum and maximum values of 0.02 and 0.4, respectively. The reduced radius at high p T allows the algorithm to reconstruct separate jets from closely spaced b-hadrons, such as in highly boosted h → bb decays.
Both the small-R and large-R jet energies are calibrated using a sequence of simulationderived corrections. Small-R jets additionally have an area-based energy subtraction applied to reduce the impact of pile-up, as well as a series of additional data-derived corrections [98].
Central small-R jets and variable-R track-jets containing b-hadrons are identified using the DL1 tagger [32]. This multivariate algorithm uses the impact parameters of ID tracks as well as information about secondary vertices and reconstructed flight paths of b-and c-hadrons within the jet. For both classes of jet a WP is chosen which tags jets containing b-hadrons with 77% efficiency in tt events. The decays of these b-hadrons can produce muons which are vetoed when building particle-flow objects and therefore not included in the energies of either the small-or large-R jets. In order to correct for this, the fourmomenta of non-isolated muons falling inside these jet cones can be added into the jet, improving the resolution of their four-momenta. For small-R (large-R) jets, this is done for the muon (two muons) closest to the jet axis. This correction is only used when calculating the mass of the Higgs boson candidate m h and was shown in ref.
[99] to improve the resolution of this measurement. Correcting m h in this way improves the ability of the fit to separate the signal from the major backgrounds.

Leptons.
Leptons are divided into 'baseline' and 'signal' categories. Events with baseline leptons are vetoed in the analysis search regions and signal leptons are used to define control regions to constrain background components. 2 The JVT selection is applied to jets with 20 GeV < p T < 60 GeV and |η| < 2.4.

JHEP11(2021)209
Electrons are reconstructed from a track which is coincident with a cluster built from energy deposits in the calorimeter [100]. They are then identified using a multivariate likelihood technique, using several features including the shape of the measured shower, the track quality and the distribution of energy within the calorimeter [101]. For this analysis, the 'LooseAndBLayer' WP [101] is used for both the signal and baseline electrons. Isolation selections are also applied to distinguish between electrons produced in the initial collision or decays of W/Z bosons or τ -leptons (prompt) and those produced in decays of other objects [101]. Requirements are placed on the energy of calorimeter clusters and the p T of tracks measured in isolation cones around the electron. For signal electrons with p T < 200 GeV and all baseline electrons the total energy of clusters within ∆R = 0.2 of the electron, excluding the electron cluster, must be less than 20% of the p T of the electron. The total p T of tracks matched to the primary vertex that lie within a cone whose size is set to the smaller of ∆R = 10 GeV/p T and 0.2, excluding the electron track, must be less than 15% of the p T of the electron. For signal electrons with p T > 200 GeV the total energy of clusters within ∆R = 0.2 of the electron is required to be less than the smaller of 0.015 × p T and 3.5 GeV. For these electrons, no track-based isolation selection is applied. Both the baseline and signal electrons are required to have |η| < 2.47. Baseline electrons are required to have p T > 7 GeV and signal electrons p T > 27 GeV. In order to ensure that they are compatible with the primary vertex, the track from which the electron is reconstructed is required to have σ(d 0 ) < 5 and |z 0 sin θ| < 0.5 mm, where σ(d 0 ) is the significance of the transverse impact parameter, z 0 is the longitudinal impact parameter, and θ is the polar angle of the track.
Muons are reconstructed by matching track segments formed in the MS to a track from the ID [102]. Identification is performed through selections on the qualities of the tracks used in the reconstruction, as well as their compatibility, for example, in the measurements of p T in the MS and ID. Similarly to electrons, isolation selections are also applied. Baseline muons are required to pass the 'Loose' identification WP and signal muons are required to pass the 'Medium' identification WP [102]. For baseline muons, the total energy of clusters within ∆R = 0.2 of the muon is required to be less than 30% of the p T of the muon. For signal muons, the total p T of tracks within ∆R = 0.2 of the muon's primary track is required to be less than 1.25 GeV. Both the baseline and signal muons are required to satisfy |η| < 2.5, σ(d 0 ) < 3 and |z 0 sin θ| < 0.5 mm. Baseline muons are required to have p T > 7 GeV and signal muons p T > 25 GeV.
Hadronically decaying τ -lepton reconstruction is seeded from R = 0.4 anti-k t jets built using the LCW-calibrated clusters [103]. As hadronic τ -lepton decays yield either one or three charged pions the jets are required to have either one or three tracks within ∆R = 0.2 of the jet axis. A recurrent neural network (RNN) classifier is used to identify the τ -leptons [33]. The inputs to the RNN are built from the clusters and tracks associated with the τ -lepton. All τ -leptons are required to pass the 'VeryLoose' WP [33] and have |η| < 2.5 and p T > 20 GeV. As there is no dedicated τ -lepton control region, no signal τ -lepton selection is defined.
Overlap removal. In order to avoid the same detector signals being interpreted as different objects, an overlap removal procedure is applied as follows. If any object is rejected -8 -

JHEP11(2021)209
at one step it is not considered in later steps. First, if any two electrons share a track the electron with the lower p T is removed. Next, any τ -leptons within ∆R = 0.2 of an electron or muon are removed. Then, any electrons which share a track with a muon are removed. If any small-R jet is within ∆R = 0.2 of an electron it is removed, and then any electron within a cone of p T -dependent size around a small-R jet is removed. If any small-R jet with fewer than three tracks has an associated muon or is within ∆R = 0.2 of one it is removed, and then any muons within a cone of p T -dependent size around a small-R jet are removed. Next, any small-R jets within ∆R = 0.2 of a τ -lepton are removed. Finally, any large-R jets within ∆R = 1.0 of an electron are removed.
Track-jets do not participate in the overlap removal as they are only used for b-tagging.
Missing transverse momentum. The missing transverse momentum (with magnitude E miss T [104]) is defined as the negative vector sum of the transverse momenta of all the observable objects in the event, plus a soft term including ID tracks matched to the primary vertex but not to any of the other objects. The E miss T reconstruction uses the baseline electrons and muons as well as all small-R jets, and employs a separate overlap removal procedure which takes into account detector signals from each object included [104]. In control regions a modified definition of E miss T is used in which electrons and muons are treated as invisible, E miss T, lep. invis. , to imitate the kinematics of the Z → νν background process. Object mismeasurements, especially of jets, are the main source of fake E miss T . Therefore, the E miss T significance (S) [35] is defined to assess the likelihood that the E miss T is really due to invisible particles or is more likely to come from mismeasurements. It is calculated using the expected resolutions of all objects which enter the E miss T calculation and the correlations between them.

Event selection
The basic target final-state topology is a Higgs boson decaying into two b-quarks produced with a significant imbalance in the measured transverse momentum. Events are divided into non-overlapping regions designed either to be enriched in the signal process (signal regions) or in a significant background process (control regions). Control regions differ from signal regions primarily through requiring the presence of one or two lepton(s), whereas signal regions veto events containing baseline leptons.
As the angle between the two b-jets produced in the Higgs boson decay is inversely proportional to the p T of the Higgs boson, in cases where the Higgs boson is significantly boosted it can become difficult to reconstruct the two b-quarks as separate jets. This motivates splitting the analysis into 'resolved' regions in which the decay products of the Higgs boson are reconstructed as two separate jets and 'merged' regions in which the entire Higgs boson decay is reconstructed as a single jet.
In b-associated production within the 2HDM+a benchmark model, the Higgs boson and DM particles are produced with an extra pair of b-quarks from gluon splitting. Therefore, to enhance sensitivity to these models, all regions are further split into those requiring exactly two b-jets and those requiring ≥3 b-jets (referred to as 2 b-tag and ≥3 b-tag, respectively).

Common selections
Events which do not have a reconstructed primary vertex are rejected. Events are also rejected if found to contain any jets with properties consistent with beam-induced backgrounds, cosmic-ray showers or noisy calorimeter cells [105].
Events are required to have E miss T > 150 GeV and are vetoed if they contain a baseline τ -lepton. In order to further reduce the background from τ -lepton decays, events are also vetoed if they have any small-R jets with ∆φ(jet, E miss T ) < 22.5 • where the track multiplicity in the jet is between 1 and 4. This selection is referred to as the 'extended τ -lepton veto'. A further source of background is E miss T arising from either leptonic heavyflavour decays in a jet or a jet which is severely mismeasured. In these cases, the E miss T tends to be aligned with the jet; therefore, events where any of the up to three leading small-R jets have ∆φ(jet, E miss T ) < 20 • are rejected. For control regions, these requirements use E miss T, lep. invis. rather than E miss T . Only loose selections are placed on the mass of the Higgs boson candidate (m h , defined in the following sections) because it is used as the discriminating variable for the final fit. The range is 50 GeV < m h < 280 GeV for the resolved regions and 50 GeV < m h < 270 GeV for the merged regions, with the lower limit chosen to be the lowest calibrated large-R jet mass and the upper limit chosen to be significantly larger than the Higgs boson mass, with the precise value being determined by the m h binning used in the fit, which depends on the available sample size.

Signal regions
The signal region selections for both the merged and resolved regions are summarised in table 1. All signal region events are required to have passed the primary E miss T trigger [44]. In order to reduce the contribution of SM processes producing E miss T through the decay W → ν, events are rejected if they contain any baseline electron or muon.

Resolved regions
The resolved regions are defined by selecting events with E miss T < 500 GeV. Events in the resolved regions are required to have at least two b-tagged small-R jets, with the two with the highest p T forming the Higgs boson candidate. The combined p T of this two-jet system (p Th ) is required to be greater than 100 GeV and its mass is corrected for nearby muons as described in section 4 to form m h .
The dominant background in the resolved region is tt production where one top quark decays leptonically, but the lepton is either not reconstructed or not correctly identified. In these cases, all the E miss T in the event (beyond that from mismeasurement) originates from the decay of one of the two W bosons, and therefore the transverse mass of the E miss In order to suppress contributions from multijet backgrounds, the object-based E miss T significance is required to satisfy S > 12. After this selection, data-driven estimates of the remaining multijet contribution to the signal regions were found to be substantially smaller than the expected statistical uncertainty of the data so the impact of multijet processes is not included in the background estimation. The studied signal models typically have fewer reconstructed jets than the dominant backgrounds. Therefore, events with exactly two b-tagged jets (2 b-tag) are required to have at most four small-R jets, where only central jets are counted. For events with at least three b-tagged jets (≥3 b-tag), this requirement is relaxed to at most five small-R jets to ensure a sufficient sample size in the corresponding control regions.

Merged regions
The merged regions are defined by selecting events with E miss T > 500 GeV. At least one large-R jet is required, and the two leading variable-R track-jets associated with the leading large-R jet are required to be b-tagged. This large-R jet is defined to be the Higgs boson candidate and its mass is corrected for nearby muons as described in section 4 to form m h . Events are separated into those which have no additional b-tagged variable-R track-jets (2 b-tag selection) and those which have at least one such b-tagged variable-R track-jet not associated with the Higgs boson candidate (≥3 b-tag selection).

Background modelling and control regions
The The tt, W +HF and Z+HF contributions are modelled using simulation with their normalisations corrected from data by using background-enriched control regions besides the signal regions. Smaller backgrounds are taken directly from simulation. These include the production of a W or Z boson in association with light jets or at most one jet containing a c-hadron, single top-quark production (dominated by production in association with a W 3 Simulated jets are labelled according to which hadrons with p T > 5 GeV are found within a cone of size ∆R = 0.3 around the jet axis. If a b-hadron is found the jet is labelled as a b-jet. If no b-hadron is found, but a c-hadron is present, then the jet is labelled as a c-jet. Otherwise the jet is labelled as a light jet. The flavour of the two leading b-tagged track-jets is used in the merged region.

JHEP11(2021)209
boson) and diboson processes. Small contributions also arise from tt processes in association with vector bosons or a Higgs boson. Another background contribution stems from vectorboson production in association with a Higgs boson (V h), which mimics the signal due to the presence of a Higgs boson peak in association with jets. In the case where the vector boson is a Z boson decaying to neutrinos this is an irreducible background. Similarly, the diboson decay ZZ → bbνν is nearly irreducible due to small difference in Z and h mass peaks (compared to the Higgs candidate mass resolution). The multijet background is negligible in all regions after the requirements on the object-based E miss T significance S and on ∆φ(jet, E miss T ), and thus not further considered. The total background estimates in all regions, including uncertainties, are determined in a simultaneous fit to all regions, which is described in section 6.
Top-quark pair production and W +HF processes contribute in the signal regions if leptons in the decays are either not identified or outside the kinematic acceptance. The main contribution arises from decays involving hadronically decaying τ -leptons. As the shape of the event variables is the same for all lepton flavours, the kinematic phase space of the signal region can be closely approximated by control regions requiring an isolated signal muon (1-muon control regions). In this case, to better approximate the signal regions which veto the presence of any leptons, E miss T, lep. invis. is used as proxy for E miss T . Also, any other variable using E miss T in its calculation, e.g. E miss T significance and m b,min/max T , is constructed using E miss T, lep. invis. . This ensures that the E miss T -related quantities in the control regions correspond to those in the signal regions. Otherwise, the 1-muon control regions are defined by the same criteria as the signal regions.
As the momentum of the Z boson does not depend on its decay mode, the main background in the signal regions, Z → νν in association with heavy-flavour jets, can be closely modelled by Z → + − events. This means that the normalisation of Z → νν+HF contribution can be corrected by measuring Z → + − +HF events. To select these events, control regions requiring exactly two baseline electrons or muons with opposite charge are defined (2-lepton control regions). These events are collected using primary triggers selecting an isolated electron or muon [106, 107]. Further, one of the electrons or muons is required to be a signal electron or muon with p T > 27 GeV or p T > 25 GeV, respectively. The invariant mass of the leptons is required to be consistent with the mass of the Z boson within 10 GeV. While keeping all other criteria of the signal regions, an additional criterion of S < 5 is imposed to suppress a remaining contribution of tt processes. To be similar to the signal regions, E miss T, lep. invis. is used as proxy for E miss T and in the calculation of any other variable using E miss T .

Statistical analysis
A binned profile likelihood fit [108, 109] is used to obtain background estimates and check the compatibility of the data with the background-only hypothesis as well as to extract upper limits at 95% confidence level (CL) on the signal cross-section. The likelihood function is constructed from a product of Poisson probability functions based on the expected signal and background yields in every region considered in the fit, as shown in  Four normalisation factors are used for the backgrounds, scaling the tt and W/Z+HF processes. Of those, two normalisation factors are used for Z boson production in association with two b-tagged jets or at least three b-tagged jets. In the 2 b-tag and ≥3 b-tag selections, tt and W +HF are normalised by a single parameter each, as the production mechanism for tt is the same in the two categories and the contribution of W +HF is minor in the ≥3 b-tag selection.
The likelihood includes all control and signal regions, as detailed in table 2, where these regions are binned to increase the sensitivity and improve the determination of the normalisation factors for tt and W/Z+HF. The signal regions are binned in the invariant mass of the Higgs boson candidate as defined in section 5. The chosen binning in the signal regions is optimised to obtain the best expected sensitivity to the signal models addressed in this paper, while also keeping statistical uncertainties in each bin low.
The 1-muon control regions are split in positive or negative muon charge in the case of the 2 b-tag regions. This improves the separation between tt and W +HF processes, as W +HF processes exhibit a charge asymmetry in pp collisions. Only inclusive event yields are used in the other regions, which are dominated by tt, and for the 2-lepton control regions.
The nominal fit results are obtained by maximising the likelihood function with respect to all parameters. Two different fit configurations are used: in the background-only profile likelihood fit, the parameters are determined in a fit to data assuming the presence of no -14 -

JHEP11(2021)209
signal. The second fit configuration instead allows for the presence of a specific signal, and is referred to as the model-dependent fit.
The test statistic q µ is constructed using the profile likelihood [110]: q µ = −2 ln(L(µ,θ µ )/L(μ,θ)), whereμ andθ are the parameters that maximise the likelihood, andθ µ are the nuisance parameter values that maximise the likelihood for a given µ. This test statistic is used to measure the compatibility of the background-only model with the observed data and to derive exclusion intervals using the CL s -prescription [111].

Systematic uncertainties
Signal and background expectations are subject to statistical, detector-related and theoretical uncertainties, which are all included in the likelihood as nuisance parameters.
Detector-related and theoretical uncertainties may affect the overall normalisation and/or shape of the simulated background and signal event distributions. Detector-related uncertainties are dominated by contributions from the jet reconstruction. Uncertainties in the jet energy scale (JES) for small-R jets [98] arise from the calibration of the scale of the jet and are derived as function of the jet p T , and also η. Further contributions emerge from the jet flavour composition and the pile-up conditions. The 'category reduction' scheme as described in ref.
[98] with 29 nuisance parameters is used. Uncertainties in the jet energy resolution (JER) depend on the jet p T and η and arise both from the method used to derive the jet resolution and from the difference between simulation and data [98], and they are included with eight nuisance parameters. Similarly, uncertainties in the jet energy resolution for large-R jets arise from the calibration, the flavour composition and the topology dependence [112]. Further uncertainties are considered for the large-R jet mass scale [112] and resolution [113].
Uncertainties due to the b-tagging efficiency for heavy-flavour jets, including c-flavour jets, are derived from tt data [114,115] and are represented by four nuisance parameters. Uncertainties are also considered for mistakenly b-tagging a light-flavour jet, with nine nuisance parameters. These are estimated using a method similar to that in ref. [116].
Uncertainties in the modelling of E miss T are evaluated by considering the uncertainties affecting the jets included in the calculation and the uncertainties in soft term's scale and resolution [117]. The pile-up in simulation is matched to the conditions in data by a reweighting factor. An uncertainty of 4% is assigned to this reweighting factor. The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [41], obtained using the LUCID-2 detector [42] for the primary luminosity measurements.
Scale factors, including their uncertainties, are calculated specifically for this analysis to correct the efficiency of E miss

JHEP11(2021)209
Modelling uncertainties impact the shape of the m h distribution, the relative acceptance between different E miss T and b-tag multiplicity bins and between signal and control regions as well as the overall normalisation of the samples that are not freely floating in the fit.
The theoretical uncertainties are dominated by modelling uncertainties in the tt and Z+HF backgrounds. For the tt and W t processes, the impact of the choice of parton shower and hadronisation model is evaluated by comparing the sample from the nominal generator set-up with a sample interfaced to Herwig 7.04 [118,119]. To assess the uncertainty in the matching of NLO matrix elements to the parton shower, the PowhegBox sample is compared with a sample of events generated with MadGraph5_aMC@NLO v2.6.2. For the W t process, the nominal sample is compared with an alternative sample generated using the diagram subtraction scheme [77,84] instead of the diagram removal scheme to estimate the uncertainty arising from the interference with tt production.
For the V +jet processes, uncertainties arising from the modelling of the parton shower and the matching scheme are evaluated by comparing the nominal samples with samples generated with MadGraph5_aMC@NLO v2.2.2. For the diboson processes, the uncertainties associated with the modelling of the parton shower, the hadronisation and the underlying event are derived using alternative samples generated with PowhegBox [60][61][62] and interfaced to Pythia 8.186 [120] or Herwig++.
For all MC samples, the uncertainties due to missing higher orders are estimated by a variation of the renormalisation and factorisation scales by a factor of two, while the PDF and α s uncertainties are calculated using the PDF4LHC prescription [121]. Table 3 gives the impact of the different sources of systematic uncertainties for selected signal models as evaluated in different model-dependent fits. The signal models with lower masses illustrate the impact of the systematic uncertainties in the resolved regions, while the models with larger mediator masses are more impacted by the merged regions. The theoretical uncertainties in the modelling of the tt background, the experimental uncertainties in the calibration of jets and the limited MC sample size show the largest impact.

Results
The post-fit background yields are determined in a background-only profile likelihood fit to data in all regions. Figure 3 shows the yields in the 1-muon and 2-lepton control regions. The post-fit normalisation factors for tt and for W +HF are found to be 0.93 ± 0.08 and 0.95 ± 0.14, respectively. For Z boson production in association with two (at least three) heavy-flavour jets the normalisation factors are determined to 1.41 ± 0.09 (1.85 ± 0.24). An upward scaling of the Z+HF background relative to the simulation was also observed in other studies [122], and was attributed to an underestimation of the g → bb rate in Sherpa. A larger scaling is observed in the region with ≥3 b-tagged jets, dominated by processes with more g → bb splittings. The uncertainty on the Z+HF normalisation factor increases in the ≥3 b-tag region due to lower statistical precision and the smaller contribution of the Z+HF background.     The distributions of the Higgs boson candidate mass m h after the background-only fit are shown in figures 4 and 5. The signal to background ratio is higher in the ≥3 b-tag signal regions, because the 2HDM+a signal model is shown with tan β = 10, where b-associated production dominates. Tables 4 and 5 present the background estimates in comparison with the observed data. No significant deviation from SM expectations is observed, with the largest deficit corresponding to a local significance of 2.3σ, and the largest excess amounting to 1.6σ. Figure 6 summarises the total yields in the signal regions as a function of E miss T . The background prediction from simulation is scaled upwards in the fit for lower E miss T values, while the simulation agrees better with the data for large E miss T values.
The results are interpreted as exclusion limits at 95% CL in the Z -2HDM and the 2HDM+a scenarios in figures 7 and 8. Considering the Z -2HDM case, Z masses up to 3 TeV are excluded for A masses of 300 GeV at 95% CL. The exclusion boundaries for the 2HDM+a scenario extend up to m a = 520 GeV for m A = 1.25 TeV for ggF production and tan β = 1. This is an improvement of about 200 GeV in m a on previous results [123], which reinterpreted the earlier h(→ bb) + E miss     Table 5. Background yields in comparison with data in the ≥3 b-tag signal regions for different E miss T ranges after a background-only fit to data. Statistical and systematic uncertainties are reported together.
The higher exclusion limit at high m A , low m a , is due to an increase of the cross-section of the a → ah process, without resonant A production. It should be noted that with the exact parameter choices adopted in this analysis, the aah coupling becomes larger than 4π for m A 1750 GeV. Moreover, as discussed in refs. [25,26] values of m A 1250 GeV (for tan β = 1) or m A 2150 GeV (for tan β = 10) would not be consistent with the requirement of having a bounded-from-below scalar potential, given the parameter choices discussed in this paper. These constraints can be relaxed substantially if the quartic couplings assume a value closer to the perturbativity limit and also in more general 2HDMs containing additional couplings as discussed in refs. [124,125]. Therefore, the above should not be considered as a strong requirement for the validity of the model predictions. At In the case of bbA production and tan β = 10, the exclusion limits extend up to m a = 240 GeV for m A = 900 GeV. The inclusion of the ≥3 b-tag tag region helps to increase the sensitivity relative to the 2 b-tag tag region by about 30-70%. The difference between observed and expected limits arises from data deficits in the ≥3 b-tag region, especially the deficit around the Higgs boson peak in the E miss T ∈ [350, 500) GeV region, as shown in figure 5. The 2HDM+a scenario with bbA production and tan β = 10 is considered for the signatures discussed in this paper for the first time.
where (A × ε) with the acceptance A and the reconstruction efficiency ε quantifies the probability for a certain event to be reconstructed within a window around the Higgs boson mass in a given signal region. The visible cross-section is obtained from the number of signal events in each signal region extracted from a fit to the m(bb) distribution as described below, divided by the integrated luminosity.
In contrast to the model-specific exclusion limits in figures 7 and 8, the upper limits on the visible cross-section are calculated without dependence on a signal model. They only assume that a resonance was produced with a mass close to 125 GeV and decays into a pair of b quarks in association with E miss T . For this purpose, the binning in m h in the signal regions is modified such that all bins under the Higgs boson peak in a range from 90 GeV to 150 GeV are merged into one bin. This bin includes the Higgs boson peak, but excludes important parts of the Z boson peak. A simultaneous fit to all control and signal regions is performed, using the modified binning. However, in order not to assume a specific signal model, a signal may only be present in the one bin under the Higgs boson peak. Although all signal regions are fitted simultaneously, the signal contributions in the different signal regions are independent of each other, which is ensured by using different signal strength parameters. leptons. The search uses proton-proton collision data recorded by the ATLAS experiment at the LHC during the data-taking periods in 2015-2018, corresponding to an integrated luminosity of 139 fb −1 . The analysis targets the 2HDM extended by dark matter particles and either a heavy vector boson Z (Z -2HDM) or a pseudoscalar singlet a (2HDM+a) as mediator particle.
Sensitivity to these models is obtained by considering different b-jet multiplicity regions, and the cases of boosted or non-boosted Higgs bosons. Further improvements relative to previous searches targeting this signature with a partial dataset from Run 2 of the LHC are obtained through better background rejection of tt processes, better identification of b-hadrons from the decay of a boosted Higgs boson, and improved suppression of multijet processes by using an object-based E miss T significance. No significant deviation from Standard Model expectations was found. Exclusion limits are set on the Z -2HDM and extend up to Z -masses of 3 TeV at 95% CL. This is an improvement of about 700 GeV relative to the analysis in ref.
[27] for heavy Z mediator masses. In the case of the 2DHM+a, the mass of the pseudoscalar a is excluded up to 520 GeV for tan β = 1 and gluon-gluon fusion production, and up to 300 GeV for tan β = 10 in b-associated production. Upper limits on the visible cross-section range from 0.05 to 3.26 fb, depending on the signal region.  [34] D. Krohn  [105] ATLAS collaboration, Selection of jets produced in 13TeV proton-proton collisions with the ATLAS detector, ATLAS-CONF-2015-029 (2015).