Search for dark matter in events with missing transverse momentum and a Higgs boson decaying into two photons in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search for dark-matter particles in events with large missing transverse momentum and a Higgs boson candidate decaying into two photons is reported. The search uses $139$ fb$^{-1}$ of proton-proton collision data collected at $\sqrt{s}=13$ TeV with the ATLAS detector at the CERN LHC between 2015 and 2018. No significant excess of events over the Standard Model predictions is observed. The results are interpreted by extracting limits on three simplified models that include either vector or pseudoscalar mediators and predict a final state with a pair of dark-matter candidates and a Higgs boson decaying into two photons.


Introduction
The discovery of a particle exhibiting the expected properties of a Standard Model (SM) Higgs boson in 2012 by the ATLAS [1] and CMS [2] Collaborations has opened up new possibilities in searches for physics beyond the SM (BSM). The precision reached since then in the production cross-section and mass measurements of the observed Higgs boson with the additional data collected by the LHC experiments [3,4] provides sharp tools to probe the possible existence of new physics.
Astrophysical data [5,6] support the existence of dark matter (DM) in our universe, while there is yet no evidence of a non-gravitational interaction between DM and SM particles, nor any indication of the microscopic nature of any possible DM particles.
Assuming the weakly interacting nature of DM, DM particles ( ) are expected to escape detection at the LHC. For this reason, searches concentrate on final states with missing transverse momentum ( miss T ) produced in association with detectable particles ( ) complementing the undetectable particles' signatures, giving rise to + miss T final states.
The detectable particle is usually chosen to be a photon [7], a or boson [8], a jet [9], a single top quark [10], or a pair of top quarks [11], all emitted from a light quark or gluon as initial-state radiation through the usual SM gauge interactions. However, can also be the neutral Higgs boson (ℎ) observed at the LHC with a mass of about 125 GeV whose radiative production is highly suppressed. In that case, collision events with ℎ produced in association with some miss T can be very sensitive probes of the structure of the BSM physics responsible for producing DM [12].
Both the ATLAS and CMS Collaborations have previously searched for such ℎ+ miss T final states using 20.3 fb −1 of collision data at √ = 8 TeV [13,14], and up to 36.1 fb −1 of collision data at √ = 13 TeV [15][16][17][18][19][20][21], considering the decay of ℎ into a pair of photons, -quarks, leptons, bosons, or bosons. In principle, the diphoton channel has lower sensitivity than the¯channel because the ℎ → branching fraction is smaller than that of ℎ →¯by two orders of magnitude. However, if a possible excess were to be seen in the -quarks final state, the diphoton channel could offer a way to cross-check this evidence for a sizeable fraction of the new-physics parameter space covered by thec hannel [22]. It also provides complementary sensitivity at lower miss T where the¯final state is limited by the experimental miss T trigger threshold, since the diphoton channel is triggered using the photon pair, allowing for much lower and better resolved miss T in the event. This paper presents an updated search for DM particles produced in association with the decay of ℎ into a pair of photons using the full LHC Run 2 collision dataset collected at √ = 13 TeV by ATLAS from 2015 to 2018. This corresponds to an integrated luminosity of 139 fb −1 , four times higher than that used in the previous analysis published by the ATLAS Collaboration [18]. Three theoretical benchmark models are considered in this analysis. The leading-order (LO) Feynman diagrams representing the production of ℎ+ miss T in these three simplified models [23] are shown in Figure 1.
In the first model, called [12], a massive vector mediator emits a Higgs boson and subsequently decays into a pair of Dirac fermionic DM candidates. DM couples to SM particles only via the new boson, and the new associated (1) baryon number symmetry ensures the stability of the DM particle in a natural way. An additional scalar particle (referred to as a baryonic Higgs boson) is introduced to break this symmetry spontaneously and generate the boson mass. The parameters of the model are: • , the boson mass; • , the coupling of the boson to the DM particle ; • , the coupling of the boson to quarks; • ℎ , the coupling between the boson and the observed Higgs boson ℎ; • sin , the mixing angle between the baryonic Higgs boson and the observed Higgs boson; and • , the mass of the fermionic dark-matter candidate .
The second and third simplified models implement different possible mediators connecting the SM spectrum to DM particles but are both derived from a general extension of the SM implementing two Higgs doublets, called 2HDM [24]. These models predict the existence of five Higgs bosons: two scalars, one being the already observed Higgs boson; one heavy pseudoscalar ; and two charged Higgs bosons ± . The first of these two models is called the -2HDM model [24], as it introduces a new vector mediator whose mass is generated by the extended Higgs sector. In this model, DM particles are produced through the decay of the pseudoscalar , giving rise to an miss T distribution that becomes harder as the mass difference between the and bosons increases. The parameters of the -2HDM model are: • , the pseudoscalar boson mass; • , the boson mass; • , the mass of the fermionic dark-matter candidate ; • tan , the ratio of the vacuum expectation values of the two Higgs doublets; • the coupling strength of the boson to quarks; and • , the mixing angle between the two neutral scalars in the 2HDM model.
The third model considered in this analysis is called 2HDM+a [25] and offers another phenomenological option where a new pseudoscalar mediator couples directly to both the SM fermions and dark-matter particles [26]. This model is of particular interest since it allows for gluon-gluon fusion production (Figures 1(c), 1(d)), which is forbidden for the mediator. This model has an extended set of 13 parameters, namely: • and , the pseudoscalar particle masses; • and ℎ , the scalar particle masses; • ± , the mass of the charged Higgs bosons; • , the fermionic DM particle mass; • , the DM Yukawa coupling; • tan , the ratio of the vacuum expectation values of the two Higgs doublets; • , the mixing angle between the two neutral scalars in the 2HDM model; • , the mixing angle between the two pseudoscalars; and • 3 , the quartic coupling of the Higgs potential, and 1 and 2 , the quartic couplings of the pseudoscalar potentials.
By choosing cos( − ) = 0 and a null coupling to quarks and leptons of the field, the 2HDM+a model can escape experimental constraints from both Higgs precision measurements and DM direct search limits.
For the particles mediating the interaction between SM and DM in these three models, it is assumed that only decays which are kinematically accessible and strictly necessary for the self-consistency of the model are considered in the decay widths [23]. In the case of the 2HDM+a model, a SM-like Higgs width is assumed for ℎ, which restricts the possible parameter space.
The analysis reported in the present paper selects events with two photons and large miss T . Searches are performed in different regions of observed diphoton transverse momentum ( T ) and event miss T significance, defined in Section 5.
The main backgrounds in the analysis correspond to either SM Higgs boson production contributions, QCD-induced non-resonant diphoton events ( , , where is a or boson), or to reducible contributions where an electron or a jet is misidentified as a photon ('fake photons') and miss T is generated either by particles escaping the detector acceptance or by neutrinos ( , +jet). An additional background contribution dominating the low miss T region is associated with resolution effects when computing the transverse energy from high-energy objects and softer contributions measured in the ATLAS calorimeters. The miss T can therefore be either 'fake', i.e. spurious values of miss T reconstructed in the detector for events with no invisible particles, or 'true', i.e. genuine miss T associated with the presence of particles escaping detection in the event.
This paper is organized as follows. Section 2 gives a brief description of the ATLAS detector. Section 3 describes the dataset and the signal and background Monte Carlo (MC) simulation samples used. Section 4 explains the reconstruction and identification of objects, while Section 5 outlines the optimization of the event selection and categorization. Section 6 summarizes the signal and background modelling. Section 7 discusses the experimental and theoretical systematic uncertainties. Section 8 presents the results and their interpretations, and finally a summary is given in Section 9.

ATLAS detector
The ATLAS detector [27][28][29] is a multipurpose particle physics detector with approximately forwardbackward symmetric cylindrical geometry. 1 The inner detector (ID) tracking system covers | | < 2.5 and consists of a silicon pixel detector, a silicon microstrip detector and a transition radiation tracker (TRT).
The ID allows precise reconstruction of charged-particle trajectories and of decay vertices of long-lived particles. The ID is surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field. A high-granularity lead/liquid-argon (LAr) sampling calorimeter measures the energy and the position of electromagnetic showers in the central (| | < 1.475) and endcap (1.375 < | | < 3.2) regions. It includes a presampler (for | | < 1.8) and three sampling layers for | | < 2.5. The longitudinal and lateral segmentation of the calorimeter allows a measurement of the shower direction without assuming that the photon originates from a specific point along the beamline. LAr sampling calorimeters with copper and tungsten absorbers are also used to measure hadronic showers in the endcap (1.5 < | | < 3.2) and forward 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upward. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the beam pipe. The pseudorapidity is defined in terms of the polar angle as = − ln[tan( /2)]. Angular distance is measured in units of Δ ≡ √︁ (Δ ) 2 + (Δ ) 2 . The photon transverse energy is T = /cosh( ), where is its energy.
(3.1 < | | < 4.9) regions, while a steel/scintillator-tile calorimeter measures hadronic showers in the central region (| | < 1.7). The muon spectrometer surrounds the calorimeters and consists of three large superconducting air-core toroid magnets, each with eight coils, a system of precision tracking chambers (| | < 2.7), and fast tracking chambers for triggering (| | < 2.4). Reconstructed events are selected by a two-level trigger system. The first-level trigger is hardware-based, while the second-level trigger is implemented in software [30].

Data and simulation samples
The collision data used in the analysis correspond to the full LHC Run 2 dataset taken by the ATLAS experiment during the period 2015 to 2018 with proton beams colliding at √ = 13 TeV. The full dataset represents an integrated luminosity of 139.0 ± 2.4 fb −1 [31] after the application of data quality requirements checking that the ATLAS detector was fully functional and that the LHC was running in stable conditions. This dataset was recorded with a mean number of about 34 interactions per bunch crossing, with a peak value of 60. Events used in this analysis were selected using a diphoton trigger requiring two reconstructed photon candidates with minimum transverse energies of 35 and 25 GeV for the leading and subleading photons, respectively, where leading (subleading) refers to the photon with the highest (second-highest) transverse energy [32]. Both photons were required by the trigger to fulfil the 'Loose' photon identification criteria in 2015 and 2016, whereas the 'Medium' criteria were used in 2017-2018 to handle the higher number of collisions per bunch crossing resulting from the higher instantaneous luminosity in those years. The trigger selections used during 2015-2018 are estimated to be fully efficient for events satisfying the offline event selection criteria presented in Section 5.
The analysis of the data sample requires the use of different types of MC generated samples to design the event selection, categorize events, and estimate systematic uncertainties applied to the statistical inference procedure used to estimate the presence of a possible DM signal in the dataset. These generated samples include the signal samples for the three models described in Section 1 and the background samples for both the non-resonant contributions ( , , ) and the SM Higgs contributions where the Higgs boson decays into a pair of photons. Due to difficulties in correctly simulating the effect of fake photons and fake miss T , the non-resonant +jet background is instead estimated using a data-driven method detailed in Section 6. The simulated non-resonant samples and data-derived +jet sample are only used to calculate uncertainties affecting the background model used in the final statistical fit.
Events from ℎ processes were modelled using the P B v2 [34-36, 41, 50] generator, which provides matrix elements at next-to-leading order (NLO) in the strong coupling constant s in the fiveflavour scheme with the NNPDF3.0nlo [51] PDF set. The functional form of the renormalization and factorization scales is set to 3 √︁ T ( ) · T (¯) · T (ℎ). 2 The generator was interfaced to P 8.2 using the A14 tune [52] and the NNPDF2.3lo [51] PDF set. The decays of bottom and charm hadrons were simulated using the E G v1.6.0 program [53].
Events from ℎ ( ℎ) processes were produced with M G 5_ MC@NLO in the four-flavour (five-flavour) scheme with the NNPDF3.0nnlo PDF [51]. The same flavour scheme is used in the matrix element calculation and the PDF. The top quark and boson decays were handled by MadSpin [54] for the correct treatment of the spin correlations of the decay products. In the case of ℎ the overlap of this process with ℎ at NLO was removed following a diagram removal technique [55,56]. The simulation of the parton shower, hadronization and underlying event was then performed by P 8.2 with the A14 tune for both the ℎ and ℎ samples.
The cross sections for the SM Higgs boson processes were calculated at NLO in electroweak theory and NNLO in QCD for the VBF, ℎ and ℎ samples [57], and next-to-next-to-next-to-leading order plus NNLL (N 3 LO+NNLL) in QCD for the ggF sample [57]. The ℎ cross section was calculated with NLO accuracy in QCD with NLO electroweak corrections [58]. The ℎ cross section was found to be negligible. All samples were normalized to the most precise available theoretical cross sections corresponding to a Higgs boson mass of 125.09 GeV [59]. The analysis assumes a branching ratio for the Higgs boson decay into two photons of 0.227% [57].
The , , and processes were simulated with the S v2.2.4 [60] generator. For the and processes, QCD NLO-accurate matrix elements for up to one parton, and LO-accurate matrix elements for up to three partons, were calculated with the Comix [61] and O L 1 [62][63][64] libraries, while for the process, QCD LO-accurate matrix elements for up to one additional parton emission were calculated. They were matched and merged with the S parton shower based on Catani-Seymour dipole factorization [61,65] using the MEPS@NLO prescription [66][67][68][69]. Samples were generated using the NNPDF3.0nnlo PDF set, along with the dedicated set of tuned parton-shower parameters developed by the S authors.
MC simulated samples for the different signal models were generated using the M G 5 [70] generator at LO accuracy, using the NNPDF3.0lo PDF set [71] for the and -2HDM signal samples or the NNPDF3.0nlo PDF set for 2HDM+a signal samples. Parton showering and hadronization were simulated using the P 8.1 [72] generator for or the P 8.2 generator for -2HDM and 2HDM+a samples, with the A14 tune and the NNPDF2.3lo PDF set [71]. Multiple samples were generated in order to scan the mediator masses and the key parameters of each model, while the default values of other parameters were set to fixed values. In general the choice of parameter values follows the recommendations of the LHC DM Forum report [23] and are based on the sensitivities expected at the LHC.
For the model, the generation scans a wide range of and DM particle masses. The mass ranges from 10 GeV to 2000 GeV while ranges from 1 GeV to 1000 GeV. The different couplings of this model are required to fulfil perturbativity bounds that allow the choice of = 1.0, = 1/3, ℎ = , and sin = 0.3 to maximize the expected cross section. The couplings choice affects only the magnitude of the cross section and not the shape of the miss T and T distributions, and in particular, a lower sin value decreases the expected yield. In this model, T and miss T grow with while the diphoton pair becomes more back-to-back with the miss T vector. Additionally, the photon pair is more collimated for higher and , since the Higgs boson is more boosted in that region of the parameter space.
2 T denotes the transverse mass of a particle, defined as T = √︃ 2 + 2 T where and T are respectively its mass and transverse momentum.
For the -2HDM model, the cross section and kinematics depend on the and masses but much less on the DM particle mass while < /2 is satisfied. Consequently, signal samples were generated for different values of ∈ [400, 1600] GeV and ∈ [200, 600] GeV and for = 100 GeV. The values of other parameters which only affect the cross section were fixed to = − /2, = 0.8, and tan = 1. In this model, miss T exhibits a peaked distribution arising from the fact that DM particles are produced from the pseudoscalar . The peak position depends strongly on the mass difference − .
For the 2HDM+a model, the parameter space is more complex than for the two first models and different scans were generated to test its richer phenomenology. This includes two sin scans corresponding to two sets of ( , ) masses, a two-dimensional scan of the -mass plane, and a two-dimensional scan of the tan -plane. The default parameter values were set to sin( − ) = 1, 3 = 1 = 2 = 3, = 1, = 10 GeV, and a decoupled Higgs boson spectrum with = = ± sizeably larger than that of the observed scalar Higgs boson with ℎ near 125 GeV.
The effects of additional collisions from the same and neighbouring bunch crossings ('pile-up') were simulated by overlaying the hard-scattering event with inelastic events generated by P 8.1 using the NNPDF2.3lo PDF set and the A3 tune [73]. Differences between the simulated and observed distributions of the number of interactions per bunch crossing were removed by applying pile-up weights to simulated events. A full simulation of the ATLAS detector [74] based on G 4 [75] was used to reproduce the detector response to SM Higgs boson processes and and backgrounds. The background and the signal samples were simulated using A F II [76], a fast simulation of the ATLAS detector response which was shown to be able to accurately simulate diphoton events.
As it is not computationally feasible to generate fully simulated samples for all interesting parameter points, an interpolation method was used to efficiently create new samples from existing fully simulated ones. First, a set of base samples covering a wide range of generator-level 'truth' values for key variables was chosen from the existing fully simulated samples. For the and -2HDM models, these variables are T and miss T , while for the 2HDM+a model, these are T , T (the transverse momentum of the system), and |Δ ( T , T )| (the difference in azimuthal angle between the and systems). For any desired additional parameter point, a new sample was generated containing only these generator-level key variables. The base samples were then reweighted on an event-by-event basis to the desired new one, using the bin-by-bin ratio of the generator-level distributions from the new sample to the corresponding distribution from the base samples. Finally, for a few validation samples, the differences between the generated and reweighted samples were used to estimate a conservative uncertainty from this procedure.

Event reconstruction
Within the ATLAS detector, photons are reconstructed from topologically connected clusters [77] of energy deposits in the EM calorimeters in the region | | < 2.37. The transition region between the barrel and endcap EM calorimeters, 1.37 < | | < 1.52, is excluded. Photon candidates matched to a conversion vertex or a track, which are consistent with originating from a photon conversion, are classified as converted photons. Those without a matched conversion vertex or track are classified as unconverted photons. The efficiency of the diphoton trigger used to select the events used in the present analysis is estimated to be greater than 99.2% on average for events passing the final event selection [32].
The calibration of the photon energy is based on a multivariate regression algorithm trained with MC samples, where the input variables are corrected with data-driven techniques. The calibrated energy is then adjusted by applying energy scale factors derived from → + − events [78]. The photon direction is reconstructed using the longitudinal segmentation of the calorimeters and a constraint from the average collision point of the proton beams. Additionally, the conversion vertex position is included in the case of converted photons.
Events are required to have at least one reconstructed primary vertex (PV), defined as a vertex associated with at least two tracks with T > 0.5 GeV. To select the correct PV of a given event, a neural network [79] that uses the pointing information from the selected photons is deployed. Although this neural-network-selected vertex is taken to be the nominally correct choice, an alternative 'hardest vertex', defined as the vertex with the highest sum of squares of the transverse momenta of associated tracks, is also used for certain miss T -related calculations described in Section 5. Photon identification is based on the lateral shower profile of the energy deposits in the first and second EM calorimeter layers and on the energy leakage fraction in the hadronic calorimeter. 'Tight' identification criteria [78] are applied, after tuning them for converted and unconverted photons separately. These criteria reduce the misidentification of hadronic jets containing large neutral components, primarily 0 particles, which decay into highly collimated photons. For transverse momenta between 30 GeV and 250 GeV, the identification efficiency for unconverted and converted photons ranges from 85% to 99%, while fake photons originating from jets have an identification efficiency of 25% to 40% [78].
To further improve fake-photon rejection, two isolation variables are defined to quantify the activity around a photon. The calorimeter-based isolation iso T is defined as the sum of the transverse energy in topological clusters of calorimeter cells within a cone of size Δ = 0.2 around the photon, correcting for the energy of the photon candidate itself as well as an average expected pile-up contribution. The track-based isolation iso T is defined as the scalar sum of the transverse momenta of all tracks with T > 1 GeV that originate from the neural-network-selected vertex and are within a cone of Δ = 0.2. Isolated photons must have iso T < 0.065 T and iso T < 0.05 T . Electrons are reconstructed from energy deposits measured in the EM calorimeter which are matched to ID tracks [78]. They are required to satisfy | | < 2.47, excluding the EM calorimeter transition region 1.37 < | | < 1.52, and have a transverse momentum T > 10 GeV. Electrons are required to satisfy the 'Medium' identification criteria based on the use of shower shape, track-cluster matching and TRT parameters [80] in a likelihood-based algorithm. Muons are reconstructed from high-quality track segments found in the muon spectrometer [81]. A matching of these segments to ID tracks is required in the region | | < 2.5. Muons are required to have | | < 2.7 and T > 10 GeV, and to satisfy the 'Medium' identification criteria [81]. Both the electrons and muons are matched to the PV via requirements on the tracks' longitudinal and transverse impact parameters, | 0 | and | 0 |. The applied requirements are | 0 | sin < 0.5 mm (where is the polar angle of the track) for electrons and muons and | 0 |/ 0 < 5(3) for electrons (muons).
Jets are reconstructed by a particle-flow algorithm [82] using noise-suppressed positive-energy topological clusters in the calorimeter [83] which are formed by the anti-algorithm [84] with a radius parameter of = 0.4. They are required to have | | < 4.4 and T > 25 GeV. To further suppress jets produced in pile-up interactions, each jet within the tracking acceptance, i.e. | | < 2.4, and with T < 60 GeV, is required to satisfy jet vertex tagger (JVT) [85] criteria used to identify the jet as originating from the PV.
To resolve ambiguities between photon, electron, muon and jet reconstruction, an overlap removal procedure is applied to avoid multiple usage of the same detector signals in the same event. The prioritization of photons in this analysis requires the removal of electrons, muons and jets within Δ = 0.4 of a selected photon. Next, jets within Δ = 0.2 of electrons are removed. In the last step, electrons and muons within Δ = 0.4 of any jet are removed.
The missing transverse momentum miss T is calculated as the magnitude of the negative vectorial sum of the transverse momenta of all selected and calibrated physics objects of an event that can be matched to the PV. A so-called 'soft term' is calculated from the residual tracks that originate from the PV but are not associated with any other object and is added to the miss T [86].

Event selection
Events are required to have at least two photon candidates within a fiducial region of the EM calorimeter defined by | | < 2.37, excluding the region of 1.37 < | | < 1.52. Photon candidates in this fiducial region are ordered according to their T and only the highest two are considered. These leading and subleading photon candidates must have T / > 0.35 and 0.25, respectively, where is the invariant mass of the two selected photons. Furthermore, events are required to have 105 < < 160 GeV. The data sideband is defined to consist of events in this region but excluding the region 120 < < 130 GeV.
In the DM production models, the Higgs boson recoils against the DM pair, resulting in large miss T in the event and large T of the diphoton system, denoted by T . The miss T distribution after the aforementioned photon requirements is shown in Figure 2. Following this preselection, a boosted decision tree (BDT) is trained using XGBoost 0. 82 [87] to discriminate between DM signals and the non-resonant diphoton background, using T and the miss T significance miss T as input variables. This latter variable is defined as the ratio of the miss T to its expected resolution: where the total transverse energy T is calculated from the scalar sum of the transverse momenta of the calibrated photons, electrons, muons, jets, and the soft term used in the miss T calculation described in Section 4. miss T is used to better control the experimental uncertainties that affect the reconstruction of the different physical objects that enter the calculation of miss T and improves the impact of miss T on the search sensitivity at low miss T . A study showed that very effective discriminating power between signal and backgrounds can be achieved by using only these two variables, considering both increases in sensitivity and in the size of systematic uncertainties associated with the use of additional variables. Figure 3 shows how different signal models considered in the present search populate the phase space of these two crucial BDT training variables, after applying both the photon requirements and the preselection. Combined with the miss T distribution in Figure 2, these figures illustrate that high-miss T backgrounds are dominated by events with true miss T from and processes, while and +jet backgrounds contribute significantly across the full diphoton candidate T range. The sum of statistical and experimental systematic uncertainties (indicated by the shaded bands) is relatively flat as a function of T but decreases in the high miss T tail, due to the large jet and miss T -related systematic uncertainties associated with the fake miss T in the and +jet components. The BDT is trained with signal taken from simulated 2HDM+a ( = 300 GeV, = 250 GeV, tan = 1, sin = 0.7) events, while the background used in training is taken from a data control region that differs from the nominal data selection by requiring that at least one photon fails either the identification or isolation requirements, with the expectation that this will not drastically change the event kinematics. The choice of training signal is motivated by its soft miss T distribution, which is relatively close to that of the background. The trained classifier performs well for all signal models, and using a different signal model with soft miss T in the training does not appreciably affect the sensitivity. For instance, for each of the four signals shown in Figure 3, the variance in sensitivity from retraining the existing classifier on a different signal is generally on the order of 10%. Figure 4 shows the distribution of the BDT score for various signals and backgrounds after training. The model trained on the data control region (dominated by +jet contributions) performs well when applied to the data sideband (dominated by true diphoton events) as the BDT is able to separate both these backgrounds from expected signals. The disagreement between the data control region and the data sideband at high BDT values is expected due to their different compositions and does not impact the final statistical results as the data control region is only used to train the BDT classifier.
Finally, events are separated into low miss T ( miss T < 150 GeV) and high miss T ( miss T > 150 GeV) regions. Such a split not only improves overall sensitivity but also complements ℎ+ miss T searches where ℎ decays into a pair of -quarks [15], which have low sensitivity in the miss T < 150 GeV region due to trigger thresholds. In each region, two categories are defined from two sequential ranges of the BDT score, with the ranges optimized to maximize the combined signal sensitivity in the two chosen categories while discarding the remaining events. The category naming scheme and corresponding ranges are summarized in Table 1, with the 'tight' categories having a higher training signal purity than the 'loose' categories.     Table 1. similar to that of the sample. Both the and templates are normalized according to their theoretical cross sections, which accounts for approximately 10% of the sideband events. The template is then normalized to approximately 70% of the sideband events, so that the sum of the , , and templates accounts for the entire true diphoton component. To represent the remaining +jet and dĳet components, a +jet template derived from a data control region in which exactly one photon fails the identification requirements is normalized to 20% of sideband events and combined with the true diphoton template.
For a given functional form, several fits are tested by varying the position of the signal peak between 121 and 129 GeV. The largest number of signal events obtained in these fits to the background-only templates is taken as Δ bkg model sig . Among the different analytic functions that are tested, the one with Δ bkg model sig smaller than 50% of the statistical uncertainty of the fitted signal yield and the least number of free parameters is chosen as the nominal background parameterization to describe the non-resonant background shape. In each category, an exponential function exp( · ) is found to be the best choice.
The non-resonant background modelling uncertainty Δ bkg model sig , which is taken as an estimate of the systematic uncertainty due to the choice of parameterization, is shown in Table 2 for each category. As a comparison of scale, its ratio to the expected non-resonant background non-res. bkg in the signal window 120 < < 130 GeV is also provided. This quantity is less than 7% in the most sensitive category, showing that the systematic uncertainty is under control in a category already dominated by statistical uncertainties. This uncertainty is implemented in the statistical model described in Section 7 as an additional signal normalized to Δ bkg model sig , times a Gaussian-constrained nuisance parameter.

Systematic uncertainties
Uncertainties from experimental and theoretical sources affect both the shapes and yields of the signal and the SM Higgs boson background, estimated from the simulated MC samples. Of these, the largest uncertainties are due to the miss T reconstruction and jets, pile-up, and signal efficiency interpolation. The theoretical systematic uncertainties include uncertainties on the factorization and renormalization scale and the parton distribution function and s (PDF+ s ). The non-resonant background is obtained directly from the fit to the data and therefore its only systematic uncertainty is the potential bias Δ bkg model sig , as described in Section 6. A summary of the experimental and theoretical uncertainties affecting the yields from SM Higgs boson processes, non-resonant background, and signal production is shown in Table 3.
The uncertainty in the integrated luminosity of the full Run 2 dataset is 1.7% [91], obtained using the LUCID-2 detector [31] for the primary luminosity measurements.
The efficiency of the diphoton trigger used to select events is evaluated in MC simulation using a trigger matching technique and in data using a bootstrap method [32]. In the diphoton invariant mass window of 105 < < 160 GeV, the trigger efficiency uncertainty affects the acceptance by 1% in each category.
The uncertainty in the vertex selection efficiency is assessed by comparing the efficiency of finding photon-pointing vertices in → + − events in data with that in MC simulation [92]. The resulting uncertainty is found to be negligible in the inclusive photon selection.
The systematic uncertainties due to the photon identification and isolation efficiencies are estimated following the prescriptions in Ref.
[78]. They are evaluated by varying the correction factors of photon selection efficiencies in MC simulation by the corresponding uncertainties. In the most sensitive category, the photon identification efficiency uncertainty is 1.3% for the SM Higgs boson and for the signal samples, while the photon isolation efficiency uncertainty is 1.4% for the SM Higgs boson and 1.3% for the signal samples. For the signal samples, an additional 2% uncertainty is added to account for photon mismodelling by the A F II simulations.
The experimental uncertainties in photon scale and resolution are obtained from Ref.
[78]. In the most sensitive category, the uncertainty in the energy scale has an effect of 1.0% on the normalization of the signals and 1.2% on the normalization of the SM Higgs boson background. The uncertainty in the energy resolution has an effect below 0.3% on the normalization of the signals and 0.4% on the normalization of the SM Higgs boson background. The effects of photon energy scale and resolution uncertainties on the signal and SM Higgs boson background mass distributions are also evaluated and parameterized in the fit to Table 3: Breakdown of the dominant systematic uncertainties. The impact of uncertainties on the yield of the SM Higgs boson processes and signal samples is shown. All production modes of the SM Higgs boson are considered together. Representative values for the impact on the most sensitive category are shown, unless one of the systematic uncertainties is not applicable to the sample, in which case the value is substituted by a "-". The "<" on the signal efficiency interpolation uncertainty indicates that the value is estimated from the maximum relative difference between the fully simulated and reweighted samples over all validation points. The impact of theoretical uncertainties is split into their effects on the migration of events between different categories and on the total cross section. Here the signal is a 2HDM+a model with = 200 GeV, = 100 GeV, tan = 1.0, sin = 0.35, which provides a conservative estimate of the size of the uncertainties for all signal points.

Source
Signals MC simulation. In particular, the impact of the scale uncertainties on the mean of these mass distributions is 0.3%, while the impact of the energy resolution on the width is 6%. An additional uncertainty in the Higgs boson mass position of 0.24 GeV from the measurement of the Higgs boson mass is considered in the fit, although this does not significantly impact the final result [59].
The migration of events among categories results from changes in object energies and momenta, due mostly to misreconstruction of jets and miss T . The experimental uncertainties in jet energy scale and resolution are propagated to the miss to the jet energy scale and resolution have an effect of 2.5% and 1.3% on the event yield of the signal and SM Higgs boson samples, respectively. The pile-up reweighting uncertainty is taken into account by propagating it through the event selection, and results in a 2.3% and 2.0% uncertainty in the event yield of the signal and SM Higgs boson samples, respectively, in the most sensitive category.
The uncertainty in the signal efficiency due to the interpolation method detailed in Section 3 is estimated by comparing the yields from fully simulated samples and reweighted samples for certain validation parameter points. The maximum relative difference between the fully simulated and reweighted samples over all these points ranges from 9% to 13%, depending on the model, and these numbers are taken as the uncertainty.
The effects of theoretical scale uncertainties on the SM Higgs boson and 2HDM+a signal samples are estimated by varying the factorization and renormalization scales up and down from their nominal values by a factor of two, recalculating the cross section in each case, and taking the largest deviation from the nominal cross section as the uncertainty. The scale uncertainties affect the event migration between the categories by 3.5% for the SM Higgs boson processes and 1.3% for the signal processes. The uncertainties in the cross section of the SM Higgs boson processes are taken from Ref. [57]. The effect of the scale uncertainty when selecting ggF Higgs boson events with high T is estimated to be 20% for each category, but the effect is negligible in the final fit because the ggF Higgs boson contribution is small. In general, the effects of theory uncertainties on the SM Higgs boson yield are small compared to the uncertainties in the non-resonant background, and do not significantly impact the sensitivity. For signals, only the uncertainties on the 2HDM+a models are taken into account, because for the and -2HDM models, no QCD vertex is calculated in the matrix element. The uncertainties in the signal cross section are not used in the fit, but are instead shown as bands on the observed limits.
For signals, PDF+ s uncertainties are estimated using the SysCalc [93] package associated with M G 5. The uncertainties in the migration of each category are estimated by varying the parameters of the NNPDF3.0lo PDF set and taking the maximum change as the uncertainty. For the SM Higgs boson, the effects of PDF+ s uncertainties on the cross sections are taken from Ref. [57]. The uncertainties in the migration of each category are estimated using the recommendations of PDF4LHC [42]. These uncertainties are 1.0% for the SM Higgs boson and 1.2% for signals. The effects of PDF+ s uncertainties on the cross section are estimated to be 2.8% for the SM Higgs boson and up to 32% for the signals. The effects of these SM Higgs boson uncertainties are once again small compared to the impact of the non-resonant background uncertainty, while the uncertainties in the signal cross sections are not used in the fit, but are instead shown as bands on the observed limits.
The uncertainty in the branching ratio ( ) of ℎ → is 1.73% [57]. The same ℎ → branching ratio is used for the SM Higgs boson and the signal models when setting limits. For the 2HDM+a models, the Higgs boson is only allowed to decay into or through a SM channel, where the ℎ → decay branching ratio is smaller than 1% when is greater than 110 GeV for = 10 GeV. This 1% effect is neglected when setting limits as it is small with respect to the uncertainty in the SM ℎ → branching ratio. Similarly, the ℎ → decay branching ratio is neglected in the model for < ℎ /2.  [95] with P for ℎ. These uncertainties are of the order of 3% for both the signal and SM Higgs boson processes.

Statistical framework
The results of the analysis are derived from a likelihood fit of the distribution in the range 105 < < 160 GeV, performed simultaneously over all four categories. The likelihood function is defined as where for each event in a category , is the observed number of events, is the expected number of events, is the value of the probability density function, are nuisance parameters, and ( ) are Gaussian constraints on the nuisance parameters.
The expected number of events is the sum of the expected yields from BSM signals, SM Higgs boson processes, the non-resonant background modelling uncertainty, and the non-resonant background: where is the signal strength, and yield and bkg model sig, represent systematic uncertainties in the resonant and non-resonant yields, as detailed in Section 7.
The nominal Higgs boson mass is set to 125.09 GeV [59] and the nominal resonant yields are fixed to values from simulation. The signal strength, non-resonant background shape parameters, and the nuisance parameters representing the systematic uncertainties are free parameters in the fit.
The test statistic is based on the likelihood ratio approach, as presented in Ref. [96].
The event yields and uncertainties in the observed data, signal, and backgrounds in the four categories within a window of 120 < < 130 GeV after the fit are shown in Table 4, for representative , -2HDM, and 2HDM+a signals. While the signal yields in each category depend on the specific model, the SM Higgs boson contribution comprises approximately 30% of the total background in the High miss T BDT tight and Low miss T BDT tight categories, due primarily to the presence of ℎ events with high genuine miss T . The uncertainties in the SM Higgs backgrounds reflect the impact of systematic uncertainties, while the uncertainties in the non-resonant background are constrained by the full diphoton mass range and the given analytic form, giving a lower uncertainty relative to the background contribution. On the other hand, the uncertainties in the signal yields reflect the fact that the fitted signal is driven mainly by the High miss T BDT tight category, where there is a large statistical uncertainty in the data. Figure 6 shows the distribution in each BDT category as well as the analytic fits to the data for the 2HDM+a model with tan = 1.0, sin = 0.35, = 600 GeV, and = 200 GeV. Different signal models have similar shapes and therefore will mainly differ in the relative signal contributions in the various categories. The fit demonstrates sensitivity to the SM Higgs boson expected production yield, with the ℎ channel providing sensitivity in the high miss T region. For each distribution, the residual between data observed in each bin and the fitted non-resonant and SM Higgs boson backgrounds is shown in the lower panel. The experimental data distribution clearly exhibits no excess with respect to the total Table 4: Event yields and uncertainties after a fit to data in the range of 120 < < 130 GeV for data, signal models, the SM Higgs boson background, and non-resonant background in each analysis category, for an integrated luminosity of 139 fb −1 . The signal samples shown correspond to a signal with = 1000 GeV and = 50 GeV, a -2HDM signal with = 800 GeV and = 500 GeV, and a 2HDM+a signal with = 600 GeV, = 200 GeV, tan = 1.0, and sin = 0. 35 background and the fit result is then interpreted to set 95% confidence level (CL) limits on the different models discussed in Section 1.

Interpretation
The observed and expected exclusion contours at 95% CL for the model in the -plane are shown in Figure 7 for sin = 0.3, = 1/3, and = 1. Compared to the results of previous ℎ+ miss T searches in the decay channel, the limit from the full LHC Run 2 dataset extends up to 1150 GeV in while it was lower than 1000 GeV in the previous ATLAS publication in this search channel [17] with early Run 2 data. The increase in the limit is only 150 GeV despite the four times larger dataset, because the signal cross section decreases strongly with increasing mass. In addition, the maximum limit on increases by more than 100 GeV to reach 280 GeV for in the range 700-950 GeV. sion, and are the numbers of protons and nucleons in the considered nucleus, set to 1 in this case. Results from the XENON [97][98][99] and DarkSide-50 [100] direct detection experiments are overlaid for the comparison. The diagonal upper branch of the limit curve reflects the fact that there are no parameters of the model that predict a cross section larger than this limit. Consequently, an observed cross section greater than this limit cannot be interpreted in this model. On the contrary, the horizontal branch delimits the sensitivity of the analysis to this model. This result improves the upper limit on the spin-independent cross section for < 2 GeV by a factor of two relative to the previous publication [17]. LHC data offers a unique window on low-mass DM candidates that complements direct DM searches in an interesting way but provides results that are more model-dependent than direct search results. The ability of each result to constrain new physics depends crucially on the model parameters. For instance, the results for the model, with model parameters sin = 0.3, = 1/3, and = 1 for this search, are more stringent than direct detection experiments for < 2 GeV and extend to DM masses well below 1 GeV, while in the case of a much lower coupling between DM and SM particles, the direct search limits may provide the more stringent constraints on this possible phenomenology. The impact of renormalization-group evolution effects [101,102] when comparing collider and direct detection limits is not considered here.
The observed and expected exclusion contours at 95% CL for the -2HDM model in the -plane are shown in Figure 9 for 0,± = , = 100 GeV, tan = 1.0, and = 0.8. The maximum limit on reaches 420 GeV for a mass of = 825 GeV. Above = 350 GeV, competing decays from →¯cause the → branching ratio to decrease quickly with increasing , resulting in the feature near = 1300 GeV. Changes in the DM mass would not significantly affect these results; such variations would only become important for very high DM masses which are not accessible to the current analysis [23]. Figures 10 and 11 show the observed and expected exclusion contours at 95% CL for the 2HDM+a model in the -and tan -planes, respectively, for sin = 0.35 and = 10 GeV. In the -scan for this benchmark point, the highest excluded is 800 GeV for = 110 GeV while the maximum excluded reaches about 260 GeV for = 600 GeV. The scan of starts from 110 GeV to avoid opening the decay of ℎ → for = 10 GeV. In the tan -scan, the region tan < 0.4 is covered with a hatched band because there the decay width of the low-mass Higgs boson is greater than 20% of its mass, which renders the cross-section calculation unreliable [22]. The shape of the limit curve closely follows the signal cross section, which is dominated by ggF for low tan and¯-fusion for high tan . Figure 12 shows the observed and expected limits on the 2HDM+a model as a function of sin for  Limits are shown at 90% CL. The results from this analysis, in which the region on the side of the hatched band inside the contour is excluded, are compared with limits from the XENON [97][98][99] and DarkSide-50 [100] experiments. The comparison is model-dependent and solely valid in the context of this model, assuming Dirac fermion DM, mixing angle sin = 0.3, and the coupling values = 1/3 and = 1. The diagonal upper branch of the limit curve reflects the fact that there are no parameters of the model that predict a cross section larger than this limit. The impact of renormalization-group evolution effects [101,102] when comparing collider and direct detection limits is not taken into consideration here. The dotted lines represent the ±1 theoretical uncertainty for the observed limit. The ±1 expected exclusion limit contour is shown as the yellow band. Above = 350 GeV, competing decays from →¯cause the → branching ratio to decrease quickly with increasing , resulting in the feature near = 1300 GeV. The dotted lines represent the ±1 theoretical uncertainty for the observed limit. The ±1 expected exclusion limit contour is shown as the yellow band. Around the threshold = 350 GeV, the competition between the resonant decay →¯and → ℎ causes the → ℎ branching ratio to decrease suddenly with a limited increase of , resulting in the feature near = 200 GeV. Above = 350 GeV, the limit is mainly driven by the increased selection efficiency due to the harder signal miss T distribution. The dotted lines represent the ±1 theoretical uncertainty for the observed limit. The ±1 expected exclusion limit contour is shown as the yellow band. The region tan < 0.4 is covered with a hatched band because there the decay width of the low-mass Higgs boson is greater than 20% of its mass, which renders the cross-section calculation unreliable [22]. The shape of the limit curve closely follows the signal cross section, which is dominated by ggF for low tan and¯-fusion for high tan .

Summary
A search for dark matter in association with a Higgs boson decaying into two photons is presented. This study is based on data collected with the ATLAS detector during 2015-2018, corresponding to an integrated luminosity of 139 fb −1 of proton-proton collisions at the LHC at a centre-of-mass energy of 13 TeV. No significant excess over the expected background is observed. Upper limits at 95% CL are set on the possible contributions to the cross sections times branching fraction of the Higgs boson decaying into two photons in association with missing transverse momentum for three models: a model, a -2HDM model, and a 2HDM+a model. The and -2HDM models provide interesting information complementary to leptophobic vector or axial-vector models, which are already largely excluded by dĳet searches [22]. Limits at 95% CL are set on the observed signal strength in the -plane for the model, the plane for the -2HDM model, and the -and tan -planes for the 2HDM+a model. A one-dimensional scan of the mixing parameter sin for , ± , = 600 GeV, = 200 GeV and tan = 1.0 in the 2HDM+a model is performed as well. Additionally, the results for the model are interpreted in terms of 90% CL limits on the DM-nucleon scattering cross section, as a function of the DM particle mass, for a spin-independent scenario. For a DM mass lower than 2 GeV, the constraint with couplings sin = 0.3, = 1/3, and = 1 placed on the DM-nucleon cross section is more stringent than limits from direct detection experiments at low DM mass, showing the complementarity between the several approaches trying to unveil the microscopic nature of DM.  The ATLAS Collaboration