Search for dark matter produced in association with a Higgs boson decaying to tau leptons at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV with the ATLAS detector

A search for dark matter produced in association with a Higgs boson in final states with two hadronically decaying τ-leptons and missing transverse momentum is presented. The analysis uses 139 fb−1 of proton-proton collision data at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV collected by the ATLAS experiment at the Large Hadron Collider between 2015 and 2018. No evidence of physics beyond the Standard Model is found. The results are interpreted in terms of a 2HDM+a model featuring two scalar Higgs doublets and a pseudoscalar singlet field. Exclusion limits on the parameters of the model in selected benchmark scenarios are derived at 95% confidence level. Model-independent limits are also set on the visible cross-section for processes beyond the Standard Model producing missing transverse momentum in association with a Higgs boson decaying into τ-leptons.


Introduction
The Standard Model (SM) of particle physics is a successful and experimentally validated theory describing elementary particles and their interactions, culminating with the discovery in 2012 [1,2] of a particle consistent with the Higgs boson (h) predicted in the SM. However, in its current form, the SM does not contain any dark matter (DM) particle candidate, whose existence is supported by astrophysical observations [3][4][5]. This provides a strong motivation for searches for physics beyond the Standard Model (BSM physics), in particular DM candidates [6].
During the Run 2 data-taking period, the Large Hadron Collider (LHC) experiments collected a substantial dataset, which can be used to search for signs of DM. There is no evidence of non-gravitational interactions between the DM and SM particles; this makes direct observation with general-purpose detectors such as ATLAS [7] or CMS [8] unlikely. The majority of the candidate DM particles χ would not interact with the material of the detector, and would escape undetected. One alternative strategy is to focus on the case of DM particles produced in association with SM particles, focusing on topologies where a single SM particle X is produced [9]. Such events do not rely on any model-specific assumptions and are characterised by a transverse-momentum imbalance (E miss T ) due to the DM particle escaping the detector, resulting in "X+E miss T " signatures, also referred to as "mono-X". The visible SM particle X can be a jet, a photon, a W or Z boson, a top quark, or a Higgs boson [10]. Typically, the couplings of such particles to light quarks and gluons of the SM are much stronger than their coupling to DM particles. The one obvious exception is a Higgs boson h whose couplings to light quarks and gluons are suppressed. This means that the mono-Higgs topologies are only significantly sensitive to scenarios in which the Higgs boson couples directly to some BSM particle participating in DM production. In such cases the direct SM-BSM interactions are probed, making the mono-Higgs searches a window into the structure of BSM physics responsible for DM production. In many theories where the Higgs couplings to BSM particles are enhanced [11], such scenarios are possible, making mono-Higgs topologies an important tool in LHC DM searches.
The One of the most straightforward and natural extensions of the SM is the two-Higgsdoublet model (2HDM) [19], extending the Higgs sector with a second scalar doublet. Such an extension of the Higgs sector is required in supersymmetric models, but is also of more general interest [20]. After electroweak symmetry breaking the model has five physical Higgs states: two neutral scalar states h and H, a pseudoscalar neutral state A, and two charged scalar states H ± . The 2HDM+a model [21] extends the baseline 2HDM by adding a new pseudoscalar singlet a with Yukawa-like coupling to the DM candidate χ and the SM fermions. The singlet acts as a mediator between the SM particles and the DM sector. This is the simplest gauge-invariant and renormalisable extension of pseudoscalar mediator DM models [22]. The model leads to h + E miss T final states through the gluon-gluon fusion and bb annihilation modes as depicted by the Feynman diagrams of figure 1.
The present analysis considers the 2HDM+a model as a benchmark and uses it to optimise the search and interpret the results. The model is described by 12 additional 1 Direct comparison with the presented search is not possible due to the difference between the models used for reinterpretation. parameters not present in the SM. Following the LHC DM Working Group's recommendations [22], two-dimensional parameter scans are performed with the following benchmark parameter choices: • Masses of the CP-odd Higgs boson A, CP-even Higgs boson H, and charged Higgs bosons H ± are set to be equal to each other.
• The mass of the fermionic DM particle χ is set to 10 GeV.
• The DM Yukawa coupling y χ to pseudoscalar a is set to 1.
• Mixing angles α and β are constrained to satisfy cos(β − α) = 0, where α is the mixing angle between CP-even Higgs bosons h and H, and tan β is the ratio of the vacuum expectation values (VEV) of the two Higgs doublets. In this so-called alignment limit the properties of the lightest Higgs boson h are very similar to those of the SM Higgs boson.
• Quartic couplings of the pseudoscalar potentials and Higgs potential λ 1P , λ 2P and λ 3 are all set to 3.
• The mass of the lightest Higgs boson h (similar to the SM Higgs boson) and the value of EW VEV v are fixed by observations to be m h ≈ 125 GeV [23] and v ≈ 246 GeV [24].
• The mixing angle between the two pseudoscalars, Θ, is required to satisfy sin Θ = 0.35 or sin Θ = 0.7.
The analysis considers two two-dimensional scans in parameter values -the first in m A and m a , the second in m A and tan β. The kinematic dependence of the signal model on m A in particular is non-trivial and the analysis is optimised for different ranges of m A ("low" and "high").
SM processes can also produce signatures with E miss T and two τ -leptons consistent with the decay products of a Higgs boson. The production of a Higgs boson in association with a Z boson decaying into neutrinos is of particular importance as, despite the relatively low cross-section, it can lead to final states indistinguishable from the signal. The production of a Z boson, a pair of vector bosons, tt, or a Higgs boson in association with a W boson can also lead to final states with E miss T and two τ -leptons. Due to the larger cross-sections, finite detector resolution and acceptance, and the ambiguity introduced by the presence of two neutrinos, these processes can mimic the potential signal events.
The paper is structured as follows. Section 2 gives a short description of the ATLAS detector. Section 3 provides information about the dataset and Monte Carlo simulations used by the analysis. Section 4 describes the various physical objects used in the analysis. In section 5, the selections applied to the dataset are summarised. Sections 6 and 7 describe the modelling and systematic uncertainties, respectively, of signal and SM backgrounds. The results of the analysis and the interpretations are discussed in section 8, with a short summary given in section 9.

ATLAS detector
The ATLAS experiment [7,25,26] is a multipurpose cylindrical particle detector at the LHC with nearly 4π coverage in solid angle. 2 At the centre of the ATLAS detector is the inner tracking detector (ID). It is surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field. The ID consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors and covers the pseudorapidity range |η| < 2.5. The next layer of the ATLAS detector is made of lead/liquidargon (LAr) sampling calorimeters that measure the energy and position of electromagnetic showers with high granularity. In the central region |η| < 1.7 the steel/scintillator-tile calorimeter is used to provide hadronic energy measurements. In the endcap and forward regions, LAr calorimeters with copper or tungsten absorbers are used for electromagnetic and hadronic shower measurements up to |η| = 4.9. The outermost layer of the ATLAS detector is the muon spectrometer. It consists of separate trigger and high-precision tracking chambers, operating in a magnetic field generated by three large air-core toroidal superconducting magnets with eight coils each.
To select and record interesting events, a two-level trigger system is used. The level-1 trigger is implemented in custom hardware. It uses information from the calorimeters and the muon spectrometer to reduce the event rate to below 100 kHz. This is followed by the high-level trigger (HLT), which is fully software-based. It is used to further reduce the event rate to 1 kHz on average.
An extensive software suite [27] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data and simulated event samples
The dataset used by the analysis was collected with the ATLAS detector during protonproton collisions at a centre-of-mass energy √ s = 13 TeV during Run 2 of the LHC (2015-2018). The recorded data correspond to an integrated luminosity of 139 fb −1 . A set of data-quality requirements are applied to the dataset to ensure that the LHC beam conditions were stable and the ATLAS detector was fully functional [28]. All data events considered by this analysis are required to pass a combined di-τ had-vis + E miss

JHEP09(2023)189
Single-top-quark processes were similarly modelled with the Powheg Box v2 [43][44][45]56] generator using the NNPDF3.0nlo set of PDFs [36]. The associated production of a top quark with a W boson (W t) and s-channel production utilised the five-flavour scheme, while the t-channel production used the four-flavour scheme. For W t processes the diagram removal scheme [93] was used to remove interference with the tt samples. The events were then interfaced to Pythia 8.230 [46], which used the A14 tune [54] and the NNPDF2.3lo set of PDFs [55]. The same configuration also applies to the production of ttH events.
The inclusive cross-sections for the previously described top-quark process samples were corrected to include higher-order effects. For tt, the cross-section was calculated to NNLO in perturbative QCD with the resummation of next-to-next-to-leading logarithmic (NNLL) soft-gluon terms using the Top++ program [47][48][49][50][51][52][53]. For W t events, the corrections were to NLO+NNLL accuracy [57,58], while s-and t-channel production was corrected to NLO accuracy [59,60]. The ttH cross-section was corrected to NLO QCD+electroweak accuracy [63].
Decays of bottom and charm hadrons for top-quark processes, including ttH, were performed by EvtGen 1.6.0 [94].
For the 2HDM+a samples, the matrix elements were generated with MadGraph5_ aMC@NLO 2.6.5 [62] at LO. Pythia 8.240 [46] with the A14 tune [54] was used to model the parton shower, hadronisation and underlying event. The CKKW-L merging procedure [95,96] was employed to match the matrix element to the parton showers. Signal Monte Carlo samples were generated with the NNPDF3.0nlo PDF set [55] and processed with a fast simulation that relies on a parameterisation of the calorimeter response [31]. Both the gluon-gluon fusion and bb annihilation signal production processes were considered. The 2HDM+a samples were generated with different choices of signal model parameters, varying two parameters at a time. In the first parameter scan, the pseudoscalar singlet mass m a was varied in the range [100,400] GeV while the mass of the pseudoscalar Higgs boson m A was varied in the range [300, 1400] GeV. The requirements tan β = 1 and sin Θ = 0.35 were applied. In the second parameter scan, m A was varied in the range [500, 1300] GeV with tan β taking values from 0.3 to 20. In this scan, constant values of m a = 250 GeV and sin Θ = 0.7 were used. The particular parameter choices are based on the recommendations from the LHC DM Working Group [22].
For all simulated samples, the effect of multiple interactions in the same and neighbouring bunch crossings (pileup) was modelled by overlaying each simulated hard-scattering event with inelastic proton-proton events generated with Pythia 8.186 [97] using the NNPDF2.3lo set of PDFs [55] and the A3 tune [98].
To compensate for differences between data and simulation, dedicated per-object correction factors derived from data are applied to jets, b-jets, τ -leptons, electrons, and muons. Additional per-event correction factors are applied to account for differences in trigger efficiency between data and simulation.

Event reconstruction
An event is considered to have a primary vertex if at least two ID tracks with p T > 500 MeV can be associated with it [99]. If there are several such vertices, the candidate with the highest sum of the squared transverse momenta of associated tracks, Σp 2 T , is considered as the hard-scattering vertex.
Jet candidates are reconstructed from particle-flow objects [100] using the anti-k t algorithm [101] with a distance parameter of R = 0.4. First, the calorimeter cells are grouped into topological clusters, with the energy measured at the electromagnetic scale to which all ATLAS calorimeters are calibrated. The energy deposited by the charged particles is then subtracted from the clusters and replaced by the momenta of matching tracks. The jet energy scale (JES) [102] calibration is applied to the jet candidates, restoring the energy to that measured at particle level. The jets are required to have a transverse momentum larger than 20 GeV and |η| < 2.5. In addition, jets with p T < 60 GeV are required to pass the tight working point of the jet-vertex-tagger (JVT) algorithm [103] to suppress background from pileup interactions.
The analysis aims for a selection orthogonal to that of the mono-Higgs search with bjets [12], so robust b-tagging with the same selection efficiency is important. Jets containing b-hadrons are identified as b-jets using the classification algorithm DL1r [104] based on deep-learning techniques. The DL1r tagger uses secondary-vertex displacement information and kinematic properties of tracks to select b-jets and reject jets originating from charm or light quarks. The b-tagging working point used for the analysis has an efficiency of 77%, as measured in an inclusive tt sample, while providing rejection factors of 4.9, 14 and 130 against c-jets, τ -jets, and light-flavour jets, respectively [105].
Electrons are reconstructed from energy deposits in the electromagnetic calorimeter that are matched to an ID track. Electron candidates are required to satisfy the loose identification criteria [106, 107], and to have p T > 10 GeV and |η| < 2.47, providing selection efficiency of 93%. Additionally, the longitudinal impact parameter z 0 has to satisfy |z 0 sin θ| < 0.5 mm to discard electron candidates not associated with the primary vertex. The objects meeting these selection criteria are called baseline electrons and participate in the computation of the missing transverse momentum and in the overlap-removal procedure, as described later. Baseline electrons that also satisfy p T > 25 GeV, fulfil the loose isolation criteria and the tight identification criteria, and satisfy the condition |d 0 |/σ(d 0 ) < 5 on the transverse impact parameter d 0 and its uncertainty σ(d 0 ) are henceforth simply referred to as electrons.
Muons are reconstructed from tracks in the muon spectrometer that are matched to tracks in the ID. Muon candidates are required to satisfy the medium identification criteria [108, 109] and to have p T > 10 GeV and |η| < 2.7. This working point has an efficiency of 97%. Additionally, the longitudinal impact parameter has to satisfy |z 0 sin θ| < 0.5 mm. The objects meeting these selection criteria are called baseline muons and participate in the computation of the missing transverse momentum and in the overlapremoval procedure. Baseline muons that also satisfy p T > 25 GeV, fulfil the loose isolation criteria, and satisfy the condition |d 0 |/σ(d 0 ) < 3 are henceforth simply referred to as muons.

JHEP09(2023)189
Reconstruction of the visible part of a hadronically decaying τ -lepton [110] is seeded by anti-k t jets with a distance parameter R = 0.4. The seed jet is built from topological clusters, calibrated with a hadronic weighting scheme [111], and is required to have p T > 10 GeV and |η| < 2.5. The τ had-vis candidate is then built from tracks and clusters within ∆R = 0.2 of the seed jet's axis. Only candidates with one or three charged tracks with a sum of ±1 are considered. τ had-vis candidates are required to have p T > 20 GeV and |η| < 2.5. Candidates in the transition region between the barrel and endcap calorimeters, 1.37 < |η| < 1.52, are excluded. The τ had-vis object's energy is calibrated with a boosted regression tree using information from the calorimeter and the particle-flow reconstruction, as well as the number of pileup interactions [112], as input. To suppress electron background, a boosted decision tree is used to reject events with electrons misidentified as τ had-vis . It is trained using tracking detector and calorimeter information as well as variables related to the ratio of the energy deposited in the calorimeter and the visible momentum of the reconstructed tracks.
To further discriminate hadronically decaying τ -leptons from background jets, a recurrent neural network (RNN) algorithm is used [113]. The algorithm uses tracking and calorimeter measurement information, as well as individual track and cluster information, as input. The analysis makes use of three different ways to define τ -lepton objects using the RNN algorithm. The first type of object (baseline τ -lepton) is required to pass the very loose identification working point. These objects are used in the overlap-removal procedure. The baseline τ -leptons that also pass the medium identification working point are simply referred to as τ -leptons. The last type of object considered by the analysis consists of τ -lepton candidates that do not pass the medium identification working point. Such objects are referred to as anti-ID τ -leptons and are used by the analysis to estimate the contribution of fake τ -leptons.
The missing transverse momentum E miss T is a measure of the transverse-momentum imbalance in the detector and is computed as the magnitude of the missing transverse momentum vector p miss T [114]. The latter is calculated as the negative vector sum of the transverse momenta of electrons, muons, τ -leptons and jets. Tracks that are associated with the primary vertex but not with any reconstructed object are also included in the calculation as the so-called "soft term".
The reconstructed objects are not exclusive, e.g. a single energy deposit in the calorimeter can be used in the reconstruction of several final-state objects. In order to resolve this ambiguity, an overlap-removal procedure is applied [115] so that at most one reconstructed object is associated with a detector signal: • If two baseline electrons share a track, only the baseline electron with the higher transverse momentum is considered.
• A baseline τ -lepton is discarded if it is closer than ∆R = 0.2 to a baseline electron or a muon.
• If a baseline muon and a baseline electron share an ID track the baseline electron is discarded. The only exception is when the baseline muon is calorimeter-tagged (i.e. a track in the ID matched to an energy deposit in the calorimeter compatible with a minimum-ionizing particle [108, 109]), then the baseline muon is discarded instead.

JHEP09(2023)189
• If a jet is closer than ∆R = 0.2 to a baseline electron or muon, it is discarded. The only exception is if the jet has three or more associated tracks, in which case it is kept instead of the baseline muon. This suppresses FSR and bremsstrahlung signatures that can be reconstructed as jets with low number of tracks.
• If any of the remaining baseline electrons or muons are closer than ∆R = 0.4 to the remaining jets, they are discarded.
• If any of the remaining jets is closer than ∆R = 0.2 to any of the remaining baseline τ -leptons, the jet is discarded.
For anti-ID τ -lepton candidates the overlap-removal procedure is performed in the same way, replacing the baseline τ -lepton identification criterion with that for anti-ID τ -leptons. However, two additional corrections are needed for the overlap-removal procedure compared to the one applied to baseline τ -leptons. First, the jets that do not pass the tight JVT working point selection requirements are given a priority over anti-ID τ -leptons not meeting the very loose identification criterion. This is done so that such jets are not masked by being treated as anti-ID τ -lepton candidates and can be vetoed successfully. The second correction is to give the same priority to b-jets. Otherwise, anti-ID τ -leptons can mask b-jets, whereas the analysis relies on an accurate counting of the number of b-jets.

Event selection
Data and simulated events are categorised in regions based on the properties of the reconstructed objects. The analysis uses three types of regions -control, validation, and signal regions. Signal regions (SRs) are regions where the sensitivity to the BSM signal model is maximised. Control regions (CRs) are regions enriched in particular background processes and are kinematically close to the SRs. A background-only fit of the simulated background processes to the data is performed in CRs using the HistFitter framework [116]. In this fit the normalisation of the dominant background processes is determined. The fitted normalisations are then propagated to the validation regions (VRs) and SRs. VRs are defined so that they lie between the SRs and CRs and serve to validate the extrapolation of background normalisation from the CRs to the SRs. The results of the background-only fit are discussed further in section 8.
Events used by the analysis are required to have a primary vertex. The events are also required to pass a jet quality criterion to reject events containing jets not originating from proton-proton collisions [117], but from beam-induced backgrounds or cosmic-ray showers.
Events are required to have exactly two τ -lepton objects that geometrically match the objects activating the di-τ had-vis part of the triggers. Both the E miss T and the di-τ had-vis parts of the trigger evolved during Run 2 to implement more efficient algorithms and adapt them to the data-taking conditions. While the HLT E miss T trigger requirement was kept constant at 50 GeV, the level-1 threshold was increased from 35 to 40 GeV in 2018. The HLT threshold for the leading-τ had-vis p T was increased from 35 GeV to 60 GeV in 2018, while the sub-leading-τ had-vis p T threshold remained constant at 25 GeV. The HLT τ had-vis -9 -

JHEP09(2023)189
Observable Selection Year 2015-2017 2018 identification algorithms were improved during Run 2, resulting in increased rejection power. The leading τ -lepton is required to have p T > 40 (65) GeV when matched to the HLT τ had-vis object passing the trigger's 35(60) GeV threshold. The sub-leading τ -lepton is required to have p T > 30 GeV. The events also have to satisfy the E miss T > 150 GeV requirement at which the trigger is operating at maximum efficiency.
Events containing an electron or a muon are vetoed. Events are required to have at most one b-jet. The common analysis preselection is summarised in table 2, defining the analysis phase space. In the following, the analysis objects in an event are assumed to be ordered by decreasing transverse momentum with indices "1" and "2" referring to the leading and sub-leading objects respectively.
In order to improve the sensitivity to the 2HDM+a signal models considered, two SRs are defined. This is done in two steps. First, a set of common SR preselection requirements is applied to suppress the dominant SM processes. After the common preselection, two non-orthogonal SRs are designed to target signal model parameter configurations with high and low masses of the heavy pseudoscalar Higgs boson A, as these yield signal events with different kinematic properties. The combination of the results from the two SRs is described in section 8.
The common SR preselection includes a veto on events with b-jets and a requirement that the two τ -leptons have opposite electric charges. The invariant mass of the two τ had-vis , or "visible invariant mass" m vis (τ 1 , τ 2 ), is constrained to be between 40 GeV and 125 GeV to ensure compatibility with a Higgs boson decay. The angular distance between the two τ -leptons is required to satisfy ∆R < 2, suppressing background events in which the two τ -leptons do not originate from a resonant decay (such as tt and W +jets). Two kinematic variables combining E miss T and the τ -leptons' p T are used to suppress Z+jets events. The event-level m tot T is defined as -10 -

JHEP09(2023)189
Common SR Preselection where p miss x and p miss y correspond to the x and y components of the missing transverse momentum vector. An event is required to have m tot T > 50 GeV to pass the common SR preselection, while m tot T < 50 GeV selection is used for trigger efficiency studies. The transverse mass of a τ -lepton is a per-object variable defined as The sum m τ 1 T + m τ 2 T is required to be larger than 100 GeV to suppress events in which E miss T is collinear with the two-τ -lepton system, as this is typical for Z boson event topologies. A summary of the common SR preselection is presented in the upper part of table 3.
A "Low m A " SR is defined on top of the common SR preselection to target signal models with m A ≤ 800 GeV. The angular distance between the two τ -leptons is tightened to be 0.6 < ∆R(τ 1 , τ 2 ) < 1.9 and the visible invariant mass of the two τ -leptons is required to be larger than 75 GeV. The transverse masses of the leading and sub-leading τ -leptons are required to be m τ 1 T > 50 GeV and m τ 2 A "High m A " SR is constructed to improve the sensitivity to models with high m A masses, leading to signal events with higher E miss T and boosted Higgs bosons. In addition to the common SR preselection, this SR requires m tot T > 400 GeV and m τ 1 T + m τ 2 T > 400 GeV. Similarly to the Low m A SR, the events in the High m A SR are separated in two m τ 1 T + m τ 2 T bins, above and below 750 GeV. -12 -

JHEP09(2023)189
The combined acceptance and efficiency of the Low m A SR bins varies from 0.0013 to 3 · 10 −6 over all of the parameter configurations of the 2HDM+a model considered by the analysis for gluon-gluon fusion production. For the High m A SR bins the combined acceptance and efficiency varies from 0.0023 to 0 because the High m A SR is only designed to be sensitive to a subset of all the parameter configurations.

Background estimation
Two sources of events are distinguished to perform the SM background estimation -those with both reconstructed τ had-vis objects originating from the hadronic decays of prompt τ -leptons (true τ -leptons), and events with at least one of the reconstructed τ had-vis objects being a non-τ -lepton object that satisfied the selection criteria, also referred to as a fake τ -lepton. Events from Z boson and multiboson production, and the majority of tt events, tend to have two true τ -leptons in the analysis phase space. Events from both W boson and multijet production, as well as approximately 25% of tt events, are expected to have at least one fake τ -lepton.
SM background processes are modelled using a combination of simulated events and data-driven methods. Events with only true τ -leptons are modelled using simulation normalised to the data in the dedicated CRs. In the following, when referring to the tt background modelling, for example, it is only the events with two true τ -leptons that are being addressed. Events with at least one fake τ -lepton are estimated in a data-driven way. Such events are referred to as fake τ -leptons or "Fake Taus" background regardless of the underlying physical process.
The normalisation of Z boson production and various processes including top-quark production (tt production, single-top-quark production, tt production in association with Z/W ) is fitted to the data in dedicated control regions. Other SM background contributions are normalised according to the theoretical predictions. For the derivation of the normalisation scale factors all top-quark-related processes are treated as one common "Top" background.
The contribution of fake τ -lepton events (i.e. events with at least one fake τ -lepton) in the control, validation and signal regions is estimated with the fake factor method [118].
The assumption is that the probability that a non-τ -lepton object (such as a jet) meets the τ -lepton identification criteria depends purely on the object's properties, such as its p T , η, number of charged-particle tracks, and origin (i.e. quark/gluon), and not on event-level variables such as E miss T or number of τ -leptons. The ratio of the number of objects that pass the medium identification working point requirement to the number of objects that do not, but instead pass the much looser anti-ID τ -lepton selection, is called the fake factor. The fake factors were measured as a function of τ -lepton's p T , |η|, number of charged-particle tracks, and origin in regions enriched with fake τ -leptons and extrapolated to the analysis phase space. In the analysis regions with the τ -lepton selection requirement replaced by that for anti-ID τ -leptons, the fake factors are applied to the anti-ID τ -leptons to obtain a prediction of the fake τ -lepton background.
-13 - The Fake Tau VR is used to validate the modelling of fake τ -leptons. The events are required to have been recorded by the di-τ + E miss T trigger, to pass the common analysis preselection as described in table 2, and to have no b-jets. To increase the fake τ -lepton purity of this VR, the two τ -leptons are required to have the same charge. The two τ -leptons are also required to have ∆R(τ 1 , τ 2 ) < 2 to remain kinematically close to the SRs. This selection ensures that the Fake Tau VR has a purity of over 85%.

JHEP09(2023)189
Two CRs are defined by inverting or relaxing some of the common SR requirements (see table 3) such that they are kinematically close to the SRs, suppress a potential signal contribution, and are enriched in Z boson and tt backgrounds respectively. Additionally, the analysis uses four VRs. Kinematically, the VRs are defined to be between the CRs and SRs and serve to validate the extrapolation of the modelling from the CRs to the SRs. The definitions of CRs and VRs are summarised in table 4 and are further discussed in the following. Figure 3 shows examples of kinematic distributions from the Top and Z CRs. The predicted background yields are shown without performing the fit and applying the resulting normalisation factors to the SM backgrounds (as described earlier in section 5). Good agreement between the predicted and observed yields is seen in all plots.
In the analysis phase space, approximately 80% of background events with a top quark arise from tt processes; the remaining events come from single-top processes. Events in the Top CR are selected using the di-τ + E miss T triggers and are required to pass the common analysis preselection summarised in table 2. To define a region enriched in tt while suppressing other common background processes such as Z+jets, the b-jet veto common to all other regions used in the analysis is relaxed and the events are required to have exactly one b-jet instead. The angular distance between the two τ -leptons is required to be ∆R(τ 1 , τ 2 ) > 1 to suppress Z boson background and to enhance the Top background purity of this CR. The two τ -leptons are also required to have opposite charges to suppress the extrapolation from the Top CR to the SRs is physically meaningful because the same underlying processes are being considered.
Events from Z boson production constitute the largest background in the analysis phase space, and a high-purity region is straightforward to define. Events in the Z CR were recorded using the same trigger as in the SRs. The events are required to pass the common analysis preselection described in table 2, and to have no b-jets. The two τ -leptons are required to have opposite charges and to satisfy ∆R(τ 1 , τ 2 ) < 2 and m vis (τ 1 , τ 2 ) < 40 GeV to improve the Z background purity. The events are also required to pass the m tot T > 50 GeV requirement. This results in the Z CR having 90% purity in Z boson background events.
While the SRs require m vis (τ 1 , τ 2 ) > 40 GeV, the Z CR demands m vis (τ 1 , τ 2 ) < 40 GeV. To validate the extrapolation to higher values of the di-τ visible invariant mass a Z VR is defined. The selection is similar to that of the Z CR, but the visible invariant mass must satisfy 40 GeV < m vis (τ 1 , τ 2 ) < 75 GeV. To ensure orthogonality to the SRs, m τ 1 T + m τ 2 T is required to be less than 100 GeV. The simulated events in the Z VR are affected by large theoretical systematic uncertainties. This is expected to be the case at large Z p T , which is the regime where this analysis operates. Systematic uncertainties related to the estimation of Z events are discussed in section 7.
The ∆R(τ 1 , τ 2 ) > 2 part of the analysis phase space has reasonably large expected Top and fake τ -lepton event yields. These two background processes are closely intertwined, considering that some of the fake τ -leptons originate from top-quark decays, and separating them is not practical. Instead, an inclusive validation region Comb VR is defined, containing Top, Fake Tau, Z boson, and diboson background processes. The Comb VR is defined by requiring that events are selected by the same trigger as used in the SRs, pass the common analysis preselection and have ∆R(τ 1 , τ 2 ) > 2. The charges of the two τ -leptons are required to be opposite and a b-jet veto is applied. The Comb VR also requires m vis (τ 1 , τ 2 ) > 125 GeV and p τ 1 +τ 2 T < 125 GeV in order to suppress potential signal contamination.
Diboson production is an important background process in the SRs. However, the overall event yields from diboson processes in the analysis phase space are too low to allow the construction of a statistically significant CR that wouldn't also be sensitive to potential signal contributions. Instead, a validation region is defined in order to validate the normalisation of the diboson backgrounds according to their cross-sections. The events in the Diboson VR are required to pass the common SR preselection described in the upper part of table 3. The events are also required to have m τ 1 T > 60 GeV, m τ 2 T > 20 GeV, and m τ 1 T + m τ 2 T > 140 GeV to increase the diboson purity of this VR. To keep the Diboson VR orthogonal to the SRs, m vis (τ 1 , τ 2 ) < 75 GeV and m tot T < 400 GeV are required. Figure 4 shows examples of kinematic distributions from the four VRs used by the analysis. The predicted background yields are shown without the derivation and extrapolation of normalisation factors in the CRs (as described in section 5). The maximum contribution from potential signal models with the parameter choices considered by the analysis does not exceed 1% in the CRs and 10% in the VRs. The only exception is the Diboson VR, where the potential signal contribution can reach 35% of the total yield.

Systematic uncertainties
The analysis considers theoretical, experimental and statistical uncertainties in signal and SM background modelling. The uncertainties are included in the likelihood fits as nuisance parameters with Gaussian probability densities. Theoretical uncertainties comprise generator-modelling-related uncertainties, cross-section uncertainties, and uncertainties related to the choice of PDF set. Such uncertainties are assumed to be uncorrelated between background processes. Systematic uncertainties in the reconstruction, identification, corrections and calibrations of various analysis objects are referred to as experimental uncertanties. These uncertainties are assumed to be correlated across processes and analysis regions.
Jet-related experimental uncertainties include uncertainties in the energy scale [119] and resolution [102], as well as jet-vertex-tagging uncertainties The uncertainty in the combined Run 2 integrated luminosity is 1.7% [123], obtained from luminosity measurements with the LUCID-2 detector [124]. The pileup reweighting applied to simulated events to match the observed conditions in data is assigned systematic uncertainties, reaching at most 6.4%. Uncertainties related to the trigger efficiency correction factors are at most 4% depending on the E miss T of the event. All of the mentioned systematic uncertainties are applied to all simulated samples.
Systematic uncertainties related to the fake-τ modelling consist of uncertainties in the number of events in the regions used to measure the fake factors and in the regions where they are applied, the statistical and systematic uncertainties of the MC samples in the regions where fake factors are applied, and the uncertainties from the extrapolation and application of fake factors. Based on a comparison of different correlation schemes the correlation effects are found to be negligible and these uncertainties are considered to be uncorrelated across regions.
Theoretical uncertainties related to variations of the PDFs are computed for all background and signal samples. The effect of using the nominal PDF and 100 variations is represented by a set of generator weights. The standard deviation of all the variations is computed and used as the PDF uncertainty. The uncertainties related to the value of the strong coupling constant α s are estimated by halving the difference of two PDF sets evaluated with α s = 0.117 and α s = 0.119. The PDF and α s systematic uncertainties are then combined in quadrature, following the PDF4LHC recommendations [75]. To estimate the uncertainties due to missing higher-order corrections, the values of renormalisation and factorisation scales µ r and µ f [125] are varied independently by a factor of 2, requiring 0.5 ≤ µ r /µ f ≤ 2. The envelope of the effects of these scale variations is used to estimate the effect of the scale uncertainty. The cross-section uncertainty is applied to all background processes that are not normalised to data in the likelihood fits.

JHEP09(2023)189
Region  Table 5. Systematic uncertainties in the post-fit SM background prediction in the six SR bins. The normalisation uncertainty includes the effect of keeping the normalisation factors as free-floating parameters in the fit. "Other" includes the uncertainties arising from trigger efficiency corrections, electrons, muons, modelling of pileup, and the E miss T computation. The individual uncertainties can be correlated and do not necessarily add in quadrature to equal the total uncertainty.
For the simulated Z and W boson samples, additional uncertainties are estimated by varying the resummation and CKKW matching scales [40,41]. For the Z(τ τ )+jets and diboson backgrounds, comparisons of the nominal Sherpa-produced samples with MadGraph5_aMC@NLO+Pythia-produced samples are used to estimate the uncertainties. For the tt and single-top-quark production processes, comparisons of the nominal Powheg+Pythia samples with MadGraph5_aMC@NLO+Pythia samples are used to estimate uncertainties related to hard scattering. Similarly, parton showering uncertainties are evaluated by comparing the Powheg+Pythia samples with Powheg+Herwig 7 samples. Uncertainties related to initial-state and final-state radiation are also considered for tt and single-top-quark production. Finally, for single-top-quark production samples, the effect of the diagram removal scheme [93] is compared with the effect of the diagram subtraction scheme [92,93] to parameterise the uncertainty in the treatment of the W t/tt interference.
For the 2HDM+a model signal samples the theoretical uncertainties related to the PDFs and renormalisation and factorisation scales are considered. Additional experimental uncertainties account for the use of the simplified parameterised calorimeter response instead of the full Geant4-based detector simulation.
A breakdown of the different sources of uncertainty in the SRs is presented in table 5. The largest contribution comes from a combination of various theoretical uncertaintiescomparisons with alternative generators for diboson, Z boson, and tt production processes, as well as PDF and scale variations. Other significant sources of uncertainty are the limited number of MC events and uncertainties in the modelling of fake τ -leptons.  The yields in each bin of the Low m A SR are summarised in table 6. Comparisons of observed and predicted yields and the corresponding statistical significances [126] are shown in figure 6 together with the benchmark parameter configurations of the signal model predictions. These signal benchmarks represent different kinematic regimes that the Low m A SR is sensitive to. The yields in the High m A SR bins are summarised in table 7 and figure 7 together with the corresponding statistical significances and benchmark parameter configurations of the signal model. No significant excess of data over the predicted SM backgrounds is observed in any of the bins. For every signal model parameter configuration considered by the analysis a signalplus-background fit is performed. The fit includes the CRs, and the Low m A or High m A SR. The contribution from the signal model is included as a free-floating parameter. The fit is performed simultaneously over all the CRs and SR bins. The probability that the observed yields in the Low m A SR or High m A SR are compatible with the signal-plusbackground hypotheses are computed. A one-sided profile-likelihood-ratio test statistic as implemented in the HistFitter framework is used. The results are obtained using the CL s prescription [127] and presented as exclusion contours at 95% confidence level (CL) in the two-dimensional parameter spaces.    The first parameter space is defined by setting tan β = 1 and sin Θ = 0.35, while varying the values of m a ∈ [100, 400] GeV and m A ∈ [300, 1200] GeV. The Low m A SR is sensitive to models with m A ≤ 800 GeV, while the High m A SR is more sensitive to the m A ≥ 1000 GeV part of the parameter space. For the Low m A SR, computations of CL s rely on the asymptotic properties of the profile-likelihood ratio [128]. For the High m A SR however, the asymptotic approximation is not valid because of the low yields in one of the regions, so Monte Carlo pseudo-experiments are used instead. The exclusion contours are derived for the Low m A SR and High m A SR separately and then combined using a method based on expected CL s values, as the two SRs are not statistically independent. The observed data yields are higher than the predicted SM backgrounds in the m τ 1 T + m τ 2 T ∈ [250, 550] GeV bins of the Low m A SR, resulting in the observed limits being weaker than the expected limits. Pseudoscalar singlets a with masses up to 300 GeV are excluded for m A = 800 GeV. The combined exclusion contour is shown in figure 8(a).

Results
The second two-dimensional parameter space considered by the analysis spans m A ∈ [500, 1300] GeV and tan β ∈ [0. 3,20] with the parameters m a and sin Θ fixed to the values of 250 GeV and 0.7 respectively. Only the Low m A SR has significant sensitivity to the signal model with these parameter value choices. Monte Carlo pseudo-experiments are used for the computation of the exclusion limits because the asymptotic approximation is not appropriate for all parameter configurations of the signal model. Since in the Low m A SR the data yields are higher than expected in the m τ 1 T + m τ 2 T ∈ [250, 550] GeV bins, the observed limits are less stringent than the expected ones for the parameter configurations of the signal model that these bins are sensitive to. This corresponds to models with lower values of m A . Signal points in the parameter space with tan β ≤ 1 are excluded for m A < 900 GeV. The exclusion contour is shown in figure 8(b).  signal-plus-background fit over all the CRs and the selected SR bin is performed. A generic signal contribution with free-floating normalisation is assumed in the SR, but not in the CRs. The profile-likelihood-ratio test statistic is evaluated using the HistFitter framework. Upper limits of 0.4−0.8 fb are placed on the visible cross-section σ vis at 95% CL. Upper limits on the number of signal events given the observed number of events and predicted number of background events, confidence levels for the background-only hypotheses, the discovery p-value, and the associated significances in all of the SR bins are summarised in table 8.

Conclusion
This paper presents a search for dark matter produced in association with a Higgs boson in final states with two hadronically decaying τ -leptons and large missing transverse momentum. The search is based on a 139 fb −1 dataset of √ s = 13 TeV proton-proton collisions recorded by the ATLAS experiment at the LHC during 2015-2018. The analysis utilises two different signal regions which in turn are divided into four and two m τ 1 T + m τ 2 T bins, respectively, to target different parts of the parameter space. No significant deviations from the Standard Model prediction is found.
The results are interpreted in the context of the 2HDM+a model where parameters are varied to set 95% CL limits in the m a -m A and m A -tan β planes. Model-independent upper limits are set on the visible cross-section σ vis for eventual beyond the Standard Model physics processes producing a large missing transverse energy in association with a Higgs -24 -

JHEP09(2023)189
boson decaying to τ -leptons. These 95% CL limits vary between 0.04 and 0.08 fb depending on which of the signal region bins is considered.
The presented work is the first exploration of the mono-Higgs signature with the Higgs boson decaying into a pair of hadronically decaying τ -leptons with the ATLAS experiment. Compared to the h + E miss T search with the Higgs boson decaying into a pair of b-quarks [12], this search offers sensitivity to the low tan β region of the m A -tan β plane. In the m a -m A plane this search has comparable sensitivity in the low m A region with an orthogonal selection, allowing for a combination of the two results.