Search for dark matter in association with a Higgs boson decaying to two photons at ﬃﬃ s p = 13 TeV with the ATLAS detector

A search for dark matter in association with a Higgs boson decaying to two photons is presented. This study is based on data collected with the ATLAS detector, corresponding to an integrated luminosity of 36 . 1 fb − 1 of proton-proton collisions at the LHC at a center-of-mass energy of 13 TeV in 2015 and 2016. No significant excess over the expected background is observed. Upper limits at 95% confidence level are set on the visible cross section for beyond the Standard Model physics processes, and the production cross section times branching fraction of the Standard Model Higgs boson decaying into two photons in association with missing transverse momentum in three different benchmark models. Limits at 95% confidence level are also set on the observed signal in two-dimensional mass planes. Additionally, the results are interpreted in terms of 90% confidence-level limits on the dark-matter – nucleon scattering cross section, as a function of the dark-matter particle mass, for a spin-independent scenario.


I. INTRODUCTION
The discovery of a particle consistent with the 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). Although strong astrophysical evidence [3,4] indicates the possible existence of dark matter (DM), there is no evidence yet for nongravitational interactions between DM and SM particles. The interaction probability of DM particles, which are produced in SM particle collisions, with a detector is expected to be tiny. Thus, many searches for DM at the Large Hadron Collider (LHC) involve missing transverse momentum (E miss T ) produced in association with detectable particles (X þ E miss T final states). In other X þ E miss T searches in proton-proton (pp) collisions, X may represent a jet or a γ=W=Z boson, which can be emitted directly from a light quark or gluon as initial-state radiation through SM gauge interactions. However, SM Higgs boson radiation from initial-state partons is highly suppressed, so events with a final state compatible with the production of a SM Higgs boson in association with E miss T can be sensitive probes of the structure of the BSM physics responsible for producing DM. The SM Higgs boson is expected to be produced from a new interaction between DM and the SM particles [5,6]. Both the ATLAS and CMS collaborations have previously searched for such topologies using considering the SM Higgs boson decay into a pair of photons or b-quarks in events with missing transverse momentum. Although the branching fraction of the SM Higgs boson decaying into a pair of photons is small, the diphoton system presented in this paper falls in a narrower mass range around the Higgs boson mass than in Ref. [11].
With the diphoton trigger, this channel is more sensitive in the low E miss T region than the channel with the SM Higgs boson decaying into a pair of b-quarks, which relies on the high E miss T trigger. This paper presents an updated search for DM particles (χ) associated with the SM Higgs boson (h) decay to a pair of photons using 36.1 fb −1 of pp collision data collected at ffiffi ffi s p ¼ 13 TeV during 2015 and 2016, where both the integrated luminosity and the center-ofmass energy are significantly higher than in the previously published ATLAS analysis [7]. Three benchmark models are considered in this analysis. The leading-order (LO) Feynman diagrams representing the production of h plus E miss T in two simplified models [12] are shown in Figs. 1(a) and 1(b). In the first model, a massive vector mediator Z 0 emits a Higgs boson and subsequently decays to a pair of Dirac fermionic DM candidates. A vector-boson mediator arises in many BSM theories through a minimal extension to the gauge sector of the SM. In scenarios where the DM couples to the SM only via the Z 0 boson [i.e., the Z 0 B model [5] represented in Fig. 1(a)], the associated U 0 ð1Þ symmetry ensures the stability of the DM particle. The baryon number B is associated with the gauge symmetry of Uð1Þ B , and an additional scalar particle (referred to as a baryonic Higgs boson) is introduced to break this symmetry spontaneously and generate the Z 0 boson mass (denoted by m Z 0 B ). The second model [from a Z 0 -two-Higgs doublet model (Z 0 -2HDM) [13], Fig. 1(b)] involves the Z 0 boson decaying to the SM Higgs boson and an intermediate heavy pseudoscalar boson A 0 , which then decays to a pair of Dirac fermionic DM particles. The minimum decay widths of the mediators are assumed for both the Z 0 B and Z 0 -2HDM models to be the sum of the partial widths for all decays into DM and quarks that are kinematically accessible [12]. Alongside those simplified models recommended in Ref. [12], a third model [referred to as the heavy-scalar model [14], Fig. 1(c)] introduces a heavy scalar boson H produced primarily via gluon-gluon fusion (ggF) with a mass in the range 2m h < m H < 2m top , in which m h and m top represent the masses of the SM Higgs boson and top quark, respectively. The upper bound on m H is introduced to avoid a large branching fraction for H → tt, which would saturate the entire H width leading to a H → hχχ branching fraction close to zero. The lower bound on m H is required to be more than twice of m h to ensure the SM Higgs boson is produced on shell. An effective quartic coupling between h, H, and χ is considered, where the DM χ is assumed to be a scalar particle. The decay branching fraction of H to h and two χ particles is assumed to be 100% for this model, to simplify the interpretations. The DM mass (m χ ) is taken to be roughly half of the SM Higgs-boson mass to ensure on-shell decay of H → hχχ, and to suppress invisible decay modes of h, as described in Ref. [15]. While no assumptions are made here as to the nature of H, it can be viewed as a part of a 2HDM þ χ scenario where H may be considered as the CP-even heavy scalar boson [14]. This paper is organized as follows. Section II gives a brief description of the ATLAS detector. Section III describes the data set and the signal and background Monte Carlo (MC) simulation samples used. Section IV explains the event reconstruction, while Sec. V outlines the optimization of the event selection and categorization. Section VI summarizes the signal and background modeling. Section VII discusses the experimental and theoretical systematic uncertainties that affect the results. Section VIII presents the results and their interpretations, and finally a summary of the results is given in Sec. IX.

II. ATLAS DETECTOR
The ATLAS detector [16,17] is a multipurpose particle physics detector with approximately forwardbackward symmetric cylindrical geometry. 1 The inner detector (ID) tracking system covers jηj < 2.5 and consists of a silicon pixel detector, a silicon microstrip detector and a transition radiation tracker (TRT). The ID allows a 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 (jηj < 1.475) and endcap (1.375 < jηj < 3.2) regions. It includes a presampler (for jηj < 1.8) and three sampling layers up to jηj < 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 beam line. LAr sampling calorimeters with copper and tungsten absorbers are also used to measure hadronic showers in the endcap (1.5 < jηj < 3.2) and forward (3.1 < jηj < 4.9) regions, while a steel/scintillator-tile calorimeter measures hadronic showers in the central region (jηj < 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 (jηj < 2.7), and fast tracking chambers for triggering (jηj < 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 [18]. ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z axis along the beam pipe. The x axis points from the IP to the center of the LHC ring, and the y axis points upward. Cylindrical coordinates ðr; ϕÞ 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Þ. The photon transverse energy is E T ¼ E= coshðηÞ, where E is its energy.

III. DATA AND SIMULATION SAMPLES
The analysis uses pp collision data with a bunch crossing interval of 25 ns, collected in 2015 and 2016 at ffiffi ffi s p ¼ 13 TeV. Only events that were recorded in stablebeam conditions, when relevant detector components were functioning properly, are considered. Events were collected using a diphoton trigger requiring two reconstructed photon candidates with transverse energies (E T ) of at least 35 and 25 GeV for the E T -ordered leading and subleading photons, respectively. The trigger efficiency with respect to the offline-reconstructed photons, measured using the same method as described in Ref. [19], is 99%. The data sample corresponds to an integrated luminosity of 36.1 fb −1 . There are on average 25 interactions in the same bunch crossing (in-time pileup) in this data sample. The MC simulation of signal and background processes is used to optimize the selection criteria, estimate uncertainties, and study the shapes of the signal and background diphoton invariant mass (m γγ ) distributions. Signal events from Z 0 B , Z 0 -2HDM, and the heavy-scalar models are generated using MADGRAPH_AMC@NLO 2.2.3 [20] at LO in quantum chromodynamics (QCD) using the NNPDF3.0LO [21] parton distribution function (PDF) set. Parton showering and hadronization are handled by the PYTHIA 8.186 [22] event generator with the A14 [23] set of tuned parameters (tune), using the NNPDF2.3LO PDF set [24]. MC samples for the Z 0 B model are generated assuming a DM mass m χ between 1 and 1000 GeV and a range of mediator masses m Z 0 B between 1 and 2000 GeV. The upper bounds of the DM and mediator masses are motivated by the limited sensitivity with the current integrated luminosity, though in principle even heavier signals can be probed with this analysis. The values of the coupling constants and mixing parameter are chosen following the recommendations of the LHC dark matter forum [12]. The couplings of the Z 0 boson to DM particles (g χ ), to quarks (g q ) and to the SM Higgs boson (g hZ 0 Z 0 =m Z 0 ) are set to 1.0, 1=3, and 1.0, respectively. The mixing angle between the baryonic Higgs boson and the SM Higgs boson is set to sin θ ¼ 0.3. The kinematic distributions predicted by the model are not strongly dependent on the coupling parameters after similar cuts either at the particle level or at the reconstruction level, and thus all samples are generated using this set of parameters. Similarly, the MC samples for the Z 0 -2HDM model are generated for ranges of values of the mediator mass m Z 0 ¼ 400 to 1400 GeV and pseudoscalar boson mass m A 0 ¼ 200 to 450 GeV for which the search is sensitive. The masses of the neutral CP-even scalar (H 0 ) and the charged scalars (H AE ) from Z 0 -2HDM model are set to 300 GeV. The DM mass m χ is set to 100 GeV, the ratio of the two-Higgs-doublet vacuum expectation values is set to tan β ¼ 1.0 and the coupling constant between the Z 0 , Higgs, and pseudoscalar bosons is set to g Z 0 ¼ 0.8, as suggested by Ref. [12]. The kinematics do not change appreciably when different values of tan β and g Z 0 are used. Moreover, in the signal process, as the DM pairs are produced as a result of the A 0 decay, there are minimal kinematic changes when the A 0 is produced onshell. For the heavy-scalar model, the events are generated with m H in steps of 10 GeV in the intervals 260 ≤ m H ≤ 270 GeV and 300 ≤ m H ≤ 350 GeV and in steps of 5 GeV in the intervals 270 ≤ m H ≤ 300 GeV for m χ ¼ 60 GeV. This mass value (m χ ¼ 60 GeV) was chosen in order to ensure on-shell decay of H → hχχ, and at the same time to suppress invisible decay modes of h. The impact of this choice on the expected sensitivity is negligible for χ mass from 50 to 70 GeV. The mass point m H ¼ 275 GeV and m χ ¼ 60 GeV is taken as an example of the kinematics that depend on the value of m H in this model.
The dominant backgrounds are resonant SM h → γγ, and nonresonant γγ, γ þ jet, Wγ, Zγ, Wγγ, and Zγγ production. For the resonant SM Higgs-boson production, events from Wh and Zh processes are generated by PYTHIA 8.186 with the A14 tune and the NNPDF2.3LO PDF set. The ggF and vector-boson fusion (VBF) samples are generated by POWHEG-BOX 2 [25][26][27][28] interfaced to PYTHIA 8.186 with the AZNLO [29] tune and the CT10 PDF set [30]. Samples of tth events are generated with MADGRAPH_AMC@NLO 2.2.3 [20] interfaced to PYTHIA 8.186 with the NNPDF3.0LO [21] PDF set. Samples of bbh events are generated by MADGRAPH_AMC@NLO 2.2.3 interfaced to PYTHIA 8.186 with the A14 tune and the NNPDF2.3LO PDF set. The nonresonant diphoton processes with associated jets are generated using SHERPA 2.1.1 [31] with the CT10 PDF set. Matrix elements are calculated with up to three partons at LO and merged with the SHERPA 2.1.1 parton shower [32] using the ME+PS@LO prescription [33]. The CT10 PDF set is used in conjunction with a dedicated parton-shower tuning developed by the SHERPA 2.1.1 authors. The Wγ, Zγ, Wγγ, Zγγ samples are generated using SHERPA 2.1.1 with the CT10 PDF set.
The normalization of nonresonant backgrounds is obtained directly from data, as described in Sec. VI. The cross sections [34] of the SM Higgs-boson processes are calculated at next-to-leading order (NLO) in electroweak theory and next-to-next-to-leading order (NNLO) in QCD for VBF, Zh and Wh samples and next-tonext-to-next-to-leading order plus next-to-next-to-leading logarithm (N 3 LO þ NNLL) in QCD for the ggF sample. The SM Higgs-boson mass is set to 125.09 GeV [35] and its branching fraction decaying into two photons is 0.227% [34]. Samples for the Z 0 B and Z 0 -2HDM models are normalized using the LO cross sections predicted by MADGRAPH_AMC@NLO 2.2.3. Samples for the heavyscalar model are normalized using the cross section of a SM Higgs boson produced in ggF at the same mass at NNLO þ NNLL in QCD.
Different pileup conditions from same and neighboring bunch crossings as a function of the instantaneous luminosity are simulated by overlaying minimum-bias events generated with PYTHIA 8.186 and EVTGEN [36] with the MSTW2008LO PDF set and the A2 [37] tune on the events of all hard processes. Differences between the simulated and observed distributions of the number of interactions per bunch crossing are removed by applying pileup weights to simulated events. Detector effects are simulated using a full simulation [38] performed using GEANT4 [39] for signals, SM Higgs processes, and Wγ, Zγ, Wγγ, and Zγγ backgrounds. The diphoton continuum background and some of the signal samples are simulated using a fast simulation based on ATLFASTII [38].

IV. EVENT RECONSTRUCTION
In each event, the observed final state is reconstructed from photons, leptons, jets, and E miss T that are built combining the related measurements provided by the various subdetectors of the experiment.
Photons are reconstructed from clusters of energy deposits in the electromagnetic calorimeter measured in projective towers. Clusters without matching tracks are classified as unconverted photon candidates. A photon is considered as a converted photon candidate if it is matched to a pair of tracks that pass a TRT-hits requirement [40] and that form a vertex in the ID which is consistent with originating from a massless particle, or if it is matched to a single track passing a TRT-hits requirement and having a first hit after the innermost layer of the pixel detector. The photon energy is calibrated using a multivariate regression algorithm trained with fully reconstructed MC samples and then corrections based on in situ techniques on data, as explained in Ref. [41]. The overall energy scale in data, as well as the difference in the constant term of the energy resolution between data and simulation, are estimated with a sample of Z boson decays to electrons recorded in 2015 and 2016. The photon direction is estimated either using electromagnetic (EM) calorimeter longitudinal segmentation (if unconverted) or conversion vertex position (if converted), together with constraints from the pp collision point.
A "tight" photon identification requirement [40] is applied to reduce the misidentification of hadronic jets containing a high-p T neutral hadron (e.g. π 0 ) decaying to two photons. The photon identification is based on the lateral profile of the energy deposits in the first and second layers of the electromagnetic calorimeter, and on the shower leakage fraction in the hadronic calorimeter. The selection requirements are tuned separately for unconverted and converted photon candidates. The identification efficiency of unconverted (converted) photons range from 85% to 95% (90% to 98%) between 25 and 200 GeV [42].
Corrections are applied to the electromagnetic showershape variables for simulated photons, to account for small differences observed between data and simulation. The diphoton mass resolution at m γγ ¼ 125 GeV is in the range 1.32-1.86 GeV [43].
To further reject hadronic backgrounds, requirements on two photon isolation variables are applied. The first variable, E iso T , is the sum of the transverse energies deposited in topological clusters [44] of cells in the calorimeter within a cone of size ΔR ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðΔηÞ 2 þ ðΔϕÞ 2 p ¼ 0.2 around each photon. The photon cluster energy and an estimate of the energy deposited by the photon outside its associated cluster are subtracted from that sum. To reduce the effects from the underlying event and pileup, the median ambient energy computed from low-p T jets in the event [45,46] is subtracted from E iso T . Good candidates are required to have an E iso T less than 6.5% of the photon transverse energy. The second variable is a track-based isolation, defined as the scalar sum of the transverse momenta of all tracks with p T > 1 GeV and consistent with originating from the primary vertex (PV) within a cone of size ΔR ¼ 0.2 around each photon. In the case of converted photon candidates, the tracks associated with the conversion are removed. Good candidates are required to have a track isolation less than 5% of the photon transverse energy. The isolation efficiency for photons, which is mostly independent of their kinematic variables, is about 90% for the ggF SM Higgs boson process.
Selected events contain at least one PV formed from two or more tracks, each with p T > 400 MeV. In each event, the diphoton PV is chosen from the PV candidates using a neural network. The input variables to this neural network are the combined z-position of the intersections of the extrapolated photon trajectories with the beam axis; the sum of the squared transverse momenta P p 2 T and the scalar sum of the transverse momenta P p T of the tracks associated with each reconstructed vertex; and the difference in azimuthal angle Δϕ between the direction defined by the vector sum of the momenta of tracks from each vertex and that of the diphoton system. Dedicated studies of Z → e þ e − are performed in order to validate the diphoton vertex identification efficiency (correct identification of the hard process vertex by the neural network) between data and simulation. The method is similar to the one used in Ref. [47]. Studies show good agreement in diphoton vertex identification efficiency between data and simulation. The efficiency to locate the diphoton PV within 0.3 mm of the production vertex is 81% for SM Higgs boson production via ggF, 82% for a heavy-scalar signal with m H ¼ 275 GeV and scalar DM m Electrons are reconstructed from energy deposits measured in the EM calorimeter which are matched to ID tracks [48]. They are required to satisfy jηj < 2.47, excluding the EM calorimeter transition region 1.37 < jηj < 1.52, and to have p T > 10 GeV. The electrons are identified using a likelihood-based algorithm that uses the track and showershape variables as input. The "medium" criteria are applied, providing an identification efficiency varying from 85% to 95% as a function of E T [49]. Loose calorimeter and track isolation requirements are applied to electrons. Pileup and the underlying event contributions in the calorimeter isolation are estimated using the ambient energy from low-p T jets and are corrected on an event-by-event basis. In the inclusive diphoton sample the combined efficiency of the isolation requirements is 98% [50].
Muons are reconstructed from high-quality track segments in the muon spectrometer. In the region jηj < 2.5, they must be matched to ID tracks. They are required to have p T > 10 GeV and jηj < 2.7 [51]. The muon "medium" criteria are applied and the identification efficiency is 96% [52]. The muon candidates must also satisfy calorimeter and track isolation criteria. The ambient energy from low-p T jets is used to correct the contributions from pileup and the underlying event. The combined isolation efficiency varies from 95% to 99% as a function of p T from 25 GeV to 60 GeV [52].
The significance of the track's transverse impact parameter with respect to the diphoton primary vertex, jd 0 j=σ d 0 , is required to be less than 5.0 for electrons and 3.0 for muons. The longitudinal impact parameter z 0 must satisfy jz 0 j sin θ < 0.5 mm for electrons and muons.
Jets are reconstructed from three-dimensional topological clusters using the anti-k t algorithm [53] with a radius parameter of R ¼ 0.4. The jets are required to have p T > 20 GeV and jηj < 4.5 for the E miss T calculation and p T > 25 GeV and jηj < 4.4 for the event selection. The jets with jηj < 2.4 and p T < 60 GeV must pass the jet vertex tagger selection [54], in which a jet is identified as originating from the diphoton primary vertex by examining the likelihood value calculated from the track information. In addition, quality criteria are applied to the jets, and events with jets consistent with noise in the calorimeter or noncollision backgrounds are vetoed [55].
The missing transverse momentum is calculated as the magnitude of the negative vectorial sum of the transverse momenta of calibrated photons, electrons, muons, and jets associated with the diphoton primary vertex. The transverse momenta of all tracks that originate from the diphoton primary vertex but are not already used in the E miss T calculation are summed and taken into account in the E miss T calculation. This term is defined as the track-based soft term [56,57]. Clusters and tracks not associated with the diphoton primary vertex are not included in the E miss T calculation, significantly suppressing the pileup effect and thus improving the E miss T resolution.

V. EVENT SELECTION
Events are required to have at least two photon candidates with p T > 25 GeV and within a fiducial region of the EM calorimeter defined by jηj < 2.37, excluding the region of 1.37 < jηj < 1.52. Photon candidates in this fiducial region are ordered according to their E T and only the first two are considered: the leading and subleading photon candidates must have E γ T =m γγ > 0.35 and 0.25, respectively, where m γγ is the invariant mass of the two selected photons. Events are required to have 105 < m γγ < 160 GeV where the diphoton mass is calculated assuming that the photons originate from the diphoton primary vertex.
In the Z 0 B and Z 0 -2HDM models of DM production, the Higgs boson recoils against the DM pair, resulting in large E miss T in the event and large p T of the diphoton candidate, denoted as p γγ T . By contrast, in the heavy-scalar model, the spectra of E miss T and p γγ T are typically shifted to smaller values. Consequently, the classification of the selected events into categories depending on E miss T , p γγ T and other kinematic quantities leads to an increase in the overall sensitivity of the analysis to these different signal models.
The background events that survive the high-E miss T requirement but are not expected to present any genuine E miss T mostly have one or several high-p T pileup jets. The misidentification of these jets therefore leads to large E miss T . These pileup jets usually originate from a pileup vertex with larger Σp 2 T than the diphoton primary vertex, where p T is the track transverse momentum associated with a single vertex. Requiring the diphoton primary vertex to be the vertex with the largest Σp 2 T in each event helps to suppress pileup effects and reject a large fraction of the fake E miss T events. Models for which ggF is the main production mode typically pass this selection, since for example the heavy scalar produced in ggF is often accompanied by radiated jets, leading to a large Σp 2 T of the vertex. An additional variable, p hard T , defined as the magnitude of the vectorial sum of the transverse momenta photons and jets in the event, provides further discrimination against events with fake E miss T in the low-E miss T region. To further reject the background events from SM Vγ and Vγγ production (where V stands for W or Z), which contribute significantly in the high-E miss T region, a lepton (electron or muon) veto is applied.
The selected events are thus divided into five categories based on the following: , where the total transverse energy P E T is calculated from the scalar sum of the transverse momenta of TABLE I. Optimized criteria used in the categorization. The categories are defined sequentially in the rows and each category excludes events in the previous row.

Category Requirements
Mono-Higgs the calibrated photons, electrons, muons and jets used in the E miss T calculation described in Sec. IV, as well as the tracks not associated with these but the PV; (ii) the diphoton transverse momentum, p γγ T ; (iii) p hard T ; (iv) the number of leptons in the event; (v) the diphoton vertex is the highest Σp 2 T vertex.
The resulting categorization scheme is shown in Table I. The categories are defined sequentially in the rows and each category excludes events in the previous row. The Z 0 B and Z 0 -2HDM signal samples are used to optimize the definition of the Mono-Higgs category, which provides most of the sensitivity to those two models. The other four categories, which provide extra sensitivity to heavy scalar boson events with softer E miss T , are optimized using simulated heavy scalar boson samples to cover the different kinematic regimes of the heavy-scalar model. Figure 2 shows the distributions of S E miss T , p hard T , and p γγ T after the selection of diphoton candidates in 120 < m γγ < 130 GeV. These distributions illustrate three kinds of signals in different S E miss T regimes, and the contributions from the different background processes. Expected distributions are shown for a Z 0 B signal with For the distributions shown in Fig. 2, the simulation is used to determine the shapes and normalizations of the Vγ and Vγγ contributions, as well as the shape of the γγ contribution, respectively. The normalizations of the γγ and γ þ jet contributions are fixed to 79% and 19% of the data yield, where these fractions are estimated from a twodimensional sideband technique by counting the number of events in which one or both photons pass or fail the identification or isolation requirements [58]. The shape of the γ þ jet contribution is determined from a data control region where one photon fails to satisfy the photon identification criteria, after subtracting the contamination expected from γγ, Vγ and Vγγ.
The slight discrepancies observed in the distributions of S E miss T , p γγ T , and p hard T in Fig. 2 do not affect the results. Discrepancies are found mainly in nonresonant backgrounds, which are estimated directly from data, as explained in Sec. VI.

VI. SIGNAL AND BACKGROUND PARAMETRIZATION
The signal and backgrounds are extracted by fitting analytic functions to the diphoton invariant mass distribution in each category. For the signal and the SM Higgsboson background, the expected normalizations are obtained from their theoretical cross sections multiplied by the product of the acceptance times efficiency from the simulation. The shapes are modeled with a double-sided crystal ball function (as defined in Ref. [43]). The shape parameters are determined by fitting the diphoton mass distribution in simulation for each category.
Both the normalization and the shape of the nonresonant background are obtained by fitting the diphoton invariant mass distribution in data for each category. A variety of analytic functions are considered for the nonresonant background parametrization: exponential functions of different-order polynomials, Bernstein polynomials of different order, and an adapted dijet function [59]. The potential bias associated with the choice of a specific analytic function to model the continuum background (referred to as the nonresonant background modeling uncertainty, ΔN nonres bkg ) is estimated for each category as the signal event yield extracted from a signal-plus-background maximum-likelihood fit to a background-only diphoton invariant mass distribution with small statistical fluctuations [43]. The backgroundonly distribution is obtained by mixing simulated γγ, γ þ jet, Vγ, and Vγγ processes. The samples of Vγ and Vγγ are weighted according to their theoretical cross section while γγ and γ þ jet samples are normalized to the number of candidates in data in the mass window 105 < m γγ < 160 GeV scaled by the fraction of each sample (79% for γγ and 19% for γ þ jet). For a given functional form, several fits are tested by varying the position of the signal peak between 115 and 135 GeV. The largest number of signal events obtained in these fits to the background-only templates is taken as ΔN nonres bkg . Among the different analytic functions that were tested, the parametrization with the smallest ΔN nonres bkg , or the minimum number of free parameters when the same ΔN nonres bkg values are obtained, is selected as the nominal background parametrization to describe the nonresonant background shape. In addition, a χ 2 test is performed to ensure that the fit is compatible with the data in each category.
The sidebands (105<m γγ <120GeV and 130<m γγ < 160GeV) of data and the simulated events are compared with each other. Due to the low statistics in the Mono-Higgs category, the potential discrepancies between data and the simulated events are obscured by the statistical uncertainties. A visible discrepancy is only observed in the high-E miss T , intermediate-E miss T , different-vertex, and rest categories with large data statistics and in these cases a data-driven reweighting is applied to the simulated events to correct the shape. In order to check whether there is a substantial improvement in χ 2 between two nested functions, several F tests [19] are performed in the data sidebands of the Mono-Higgs category where the simple exponential function is chosen. The F tests are performed by comparing the simple exponential function to other higher-order functions. A test statistic F is calculated by using the resulting values of χ 2 , and its probability is compared with that expected from a Fisher distribution with the corresponding number of degrees of freedom. The hypothesis that a higher-order function is not needed is rejected if the F-test probability is less than 5%. The tests show that the higher-order functions do not provide a significantly better fit and the lower-order function is sufficient. The selected analytic function along with its nonresonant background modeling uncertainty in each category, which is taken as an estimate of the systematic uncertainty due to the choice of parametrization, is shown in Table II.

VII. SYSTEMATIC UNCERTAINTIES
Uncertainties from experimental and theoretical sources affect the signal efficiency and the SM Higgs-boson background yield estimated from the simulated MC samples. The nonresonant background is obtained directly from the fit to the data. The only systematic uncertainty in the nonresonant background is the potential bias ΔN nonres bkg . A summary of the experimental and theoretical uncertainties with respect to the yield of the background from SM Higgs-boson processes, nonresonant background, and signal production is shown in Table III.
The uncertainty in the combined 2015 and 2016 integrated luminosity is 3.2%. It is derived, following a methodology similar to that detailed in Ref.
[60], from a calibration of the luminosity scale using x-y beam-separation scans performed in August 2015 and May 2016.
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 [18]. In the diphoton invariant mass window of 105 < m γγ < 160 GeV, the trigger efficiency uncertainty is found to be 0.4%.
The uncertainty in the vertex selection efficiency is assessed by comparing the efficiency of finding photonpointing vertices in Z → e þ e − events in data and MC simulation [19], for which each electron track is removed TABLE II. The analytic functions used to model the nonresonant m γγ distribution, the uncertainty in the signal due to the nonresonant background modeling (ΔN nonres bkg ) and the relative uncertainties in the expected nonresonant background and signal per category. The relative uncertainty in the signal is evaluated in the Mono-Higgs category by using the Z 0 B yields and in the other categories by using the heavy-scalar yields as N signal from  28   TABLE III. Breakdown of the dominant systematic uncertainties in the range of 105 < m γγ < 160 GeV. The uncertainties in the yield of signals, the background from SM Higgs-boson processes, and nonresonant background are shown. All production modes of the SM Higgs boson are considered together. Values for the impact on all categories are shown, unless one of the systematic uncertainties is not applicable to the sample, in which case the value is substituted by a "Á Á Á". If a given source has a different impact on the various categories, the given range corresponds to the smallest and largest impacts among categories or among the different signal models used in the analysis. In addition, the potential bias coming from nonresonant background mismodeling is shown relative to the background. and its cluster is recalibrated as a photon cluster. The efficiency of this selection in data is found to be in agreement with the simulation within 0.01%. The systematic uncertainties due to the photon energy scale and resolution are adapted from run-1 results [41], with minor updates in case of data-driven corrections using the run-2 data. The uncertainty in the energy scale is less than 1% in the p T range of the photons used in this analysis; the uncertainty in the energy resolution is smaller than 2%.
Uncertainties in photon identification and isolation efficiencies are estimated, and their impact on the number of events in each category is quantified. The uncertainty in the photon identification efficiency [42] is calculated from the difference between applying and not applying the corrections to the electromagnetic-shower-shape variables in simulation. The resulting uncertainty in the photon identification efficiency is lower than 3.8% for SM Higgs background in all categories, 2.9% for simplified model samples, and 4.3% for the heavy-scalar model. The uncertainty in the photon calorimeter isolation efficiency is calculated from efficiency differences between applying and not applying corrections derived from inclusive photon events to the isolation variables in simulation. The measurements of the efficiency correction factors using inclusive photon events are used to derive the efficiency uncertainty in the photon track isolation uncertainty. The photon isolation efficiency uncertainty is found to be smaller than 1.6% for the SM Higgs background and 1.2% for all signal samples.
Migration of events among categories occurs owing to changes in the energy of identified particles and jets, mostly due to the misreconstruction of jets and E miss T . The uncertainties in jet energy scale, resolution, and jet vertex tagger are propagated to the E miss T calculation. In addition, the uncertainties in the scale and resolution of the E miss T soft term are estimated in 2016 data using the method described in Ref. [61]. The overall jet and E miss T uncertainties in the SM Higgs-boson processes are 6, 8, 23, 22, and 1% for the Mono-Higgs, high-E miss T , intermediate-E miss T , differentvertex, and rest categories, respectively. For signal processes, the overall jet and E miss T uncertainties range from 1.0% to 1.4%. The choice of the diphoton vertex also affects the E miss T reconstruction. It introduces an additional uncertainty derived from the data-to-MC comparison in Z → e þ e − events. This systematic uncertainty affects the processes with no genuine E miss T and is estimated in each category. For the SM Higgs-boson production, it is found to be 0.5% in the Mono-Higgs category and up to 1.9% in the other categories. It is less than 0.1% for signal processes. Requiring the diphoton primary vertex to be the vertex with the largest Σp 2 T in the event introduces an uncertainty of about 6% in the SM Higgs-boson production in the high-E miss T and the intermediate-E miss T categories and 1.8% in the heavy-scalar signals in those categories. This systematic uncertainty differs between the signal and background because there is a correlation between the E miss T and the pileup vertex Σp 2 T for different processes. For other signal samples, these uncertainties are at most a few percent in any category. The pileup reweighting uncertainty is taken into account by propagating it through the event selection, and results in a 0.2%-5.6% uncertainty in the event yield of the signal and SM Higgs-boson samples.
The nonresonant background contribution is not affected by the same systematic uncertainties that characterize the MC simulation, since it is extracted from the data. The potential bias is added as a systematic uncertainty to account for the possible impact on the fitted signal yields of nonresonant background mismodeling. It is the leading systematic uncertainty of the analysis but is only one-third as large as the statistical uncertainty. The ratio of the potential bias in the signal yield to the expected nonresonant background (ΔN nonres bkg =N nonres bkg ) in the range 120 < m γγ < 130 GeV is 9.8% in the Mono-Higgs category, 4.0% in the high-E miss T category, 1.3% in the intermediate-E miss T category, 0.5% in the different-vertex category, and 0.1% in the rest category.
The predicted cross sections of the SM Higgs-boson and signal processes are affected by uncertainties due to missing higher-order terms in perturbative QCD. These uncertainties are estimated by varying the factorization and renormalization scales up and down from their nominal values by a factor of 2, recalculating the cross section in each case, and taking the largest deviation from the nominal cross section as the uncertainty. The uncertainty related to the renormalization and factorization scales is 0.6%-11% for signal and 2.5%-6.0% for the SM Higgs-boson processes [34]. For the three signal processes, the effect of PDF þ α S uncertainties including acceptance and normalization is 11%-25%. It is estimated using the recommendations of PDF4LHC [34]. Both intra-PDF and inter-PDF uncertainties are extracted. Intra-PDF uncertainties are obtained by varying the parameters of the NNPDF3.0LO PDF set, while inter-PDF uncertainties are estimated using alternative PDF sets (CT14 [62] at LO and MMHT2014 [63] at LO). The final inter-PDF uncertainty is the maximum deviation among all the variations from the central value obtained using the NNPDF3.0LO PDF set. In the case of the SM Higgs-boson processes, the effect of α S and the choice of PDFs range from 2% to 6%, which are taken from Ref. [34]. The uncertainty in the branching fraction (B) of h → γγ is 1.73% [34]. The uncertainty in the multiple parton-parton interactions is estimated by switching them on and off in PYTHIA 8.186 in the production of the ggF SM Higgs-boson sample. The resulting uncertainty in the number of events in this sample conservatively reaches 0.4% in the rest category, 5.8% in the differentvertex category, 3.8% intermediate-E miss T and differentvertex category, 3.4% in the high-E miss T category, and 1.4% in the Mono-Higgs category.

VIII. RESULTS
The results for the analysis are derived from a likelihood fit of the m γγ distribution in the range 105 < m γγ < 160 GeV. The SM Higgs boson mass is set to 125.09 GeV [35]. The impact due to the SM Higgs-boson mass uncertainty is negligible. The signal strength, the background shape parameters, and the nuisance parameters representing the systematic uncertainties are set to be free parameters. The SM Higgs yields are taken from the SM predictions. The systematic uncertainty of each nuisance parameter is taken into account by multiplying the likelihood by a Gaussian penalty function centered on the nominal value of this parameter with a width set to its uncertainty. The nominal value of each SM Higgs-boson background nuisance parameter (including its yield) is taken from the simulation normalized to the SM theoretical predictions. Since adding all the other categories to the Mono-Higgs category only brings a 1% sensitivity improvement for both the Z 0 B and Z 0 -2HDM signals, the results are only obtained from this category for these two simplified models. In contrast, results for the heavy-scalar model are obtained from a simultaneous fit of all the categories.
The event yields in the observed data, expected signal and backgrounds in the five categories within a window of 120 < m γγ < 130 GeV are shown in Table IV. The signal samples shown correspond to a Z 0 B signal with m Z 0 B ¼ 200 GeV and m χ ¼ 1 GeV, a Z 0 -2HDM signal with m Z 0 ¼ 1000 GeV, m A 0 ¼ 200 GeV, and m χ ¼ 100 GeV, and a heavy-scalar signal with m H ¼ 275 GeV and m χ ¼ 60 GeV. For each benchmark signal model, the selection efficiency times acceptance denoted by "A × ϵ" is also shown. The yields for the nonresonant background are obtained from a fit to data while SM Higgs-boson events are estimated from the simulation. The uncertainties correspond to the quadrature sum of the statistical and systematic uncertainties.

A. Limits on the visible cross section
The observed yields agree with the SM background predictions, as shown in Table IV, and no significant excess of events is observed. Upper limits are set on the visible cross section σ BSM vis ≡ ðA × ϵ × σ × BÞ BSM for BSM physics processes producing missing transverse momentum and a SM Higgs boson decaying into two photons. The SM background prediction is excluded from this BSM visible cross section. Table V shows the observed and expected 95% confidence level (C.L.) upper limits on σ BSM vis , which are calculated using a one-sided profile-likelihood ratio and the C:L: s formalism [64,65] with the asymptotic approximation in Ref. [66]. The same parametrizations for the BSM signal and the total SM Higgs-boson background are used in each of the five different categories. The fits are performed individually in each category. The statistical uncertainty is dominant. The systematic uncertainties worsen the limits by about 10% (7% from the nonresonant TABLE IV. Event yields in the range of 120 < m γγ < 130 GeV for data, signal models, the SM Higgs-boson background and nonresonant background in each analysis category, for an integrated luminosity of 36.1 fb −1 . The signal samples shown correspond to a heavy-scalar sample with m H ¼ 275 GeV and m χ ¼ 60 GeV, a Z 0 B signal with m Z 0 B ¼ 200 GeV and m χ ¼ 1 GeV and a Z 0 -2HDM signal with m Z 0 ¼ 1000 GeV, m A 0 ¼ 200 GeV, m H 0;AE ¼ 300 GeV and m χ ¼ 100 GeV. For each benchmark signal model, the selection efficiency times acceptance denoted as "A × ϵ" is also shown. The yields for nonresonant background are obtained from a fit to data while SM Higgs-boson events are estimated from the simulation. The uncertainties correspond to the quadrature sum of the statistical and systematic uncertainties.   B. Interpretations of the Z 0 B and Z 0 -2HDM models Figure 3 shows the m γγ distributions in the Mono-Higgs category as well as the fits for a Z 0 B benchmark point with m Z 0 B ¼ 200 GeV and m χ ¼ 1 GeV. No significant excess is observed in this category. Upper limits are set on the production cross sections in the two theoretical models considered. Figure 4(a) shows the observed and median expected 95% C.L. upper limits on σðpp → hχχÞ × Bðh → γγÞ as a function of the mediator mass m Z 0 B for a DM mass of 1 GeV. The cross sections times branching fraction of h → γγ larger than 2.3 fb are excluded for the full range of m Z 0 B between 10 and 2000 GeV at 95% C.L., and the Z 0 B model is excluded with Z 0 B masses below 850 GeV for a DM mass of 1 GeV.
In the Z 0 -2HDM scenario, the observed and median expected 95% C.L. upper limits on σðpp → hχχÞ× Bðh → γγÞ are shown in Fig. 4(b), as a function of the pseudoscalar boson mass m A 0 for m χ ¼ 100 GeV and m Z 0 ¼ 1000 GeV. The masses of the neutral CP-even scalar (H 0 ) and the charged scalars (H AE ) from Z 0 -2HDM model are set to 300 GeV. The theoretical cross section starts from m A 0 ¼ 201 GeV. The working point with m A 0 ¼ 200 GeV is excluded since the resonant production of DM particles at m χ ¼ 100 GeV significantly increases the cross section of the process. To avoid the resonant regime where m A 0 ¼ 200 GeV and m χ ¼ 100 GeV and allow a better limit interpolation, the point m A 0 ¼ 201 GeV is shown in this plot instead of 200 GeV. The drop of the theoretical prediction at m A 0 ¼ 345 GeV is due to a rapid change in the width when A 0 decaying to tt is kinematically TABLE V. Observed and expected upper limits (at 95% C.L.) on the visible cross section for BSM physics processes producing missing transverse momentum and a SM Higgs boson decaying into two photons. Limits are presented for the five different categories. The AE1σ exclusion from the expected limits are also given. For all the simulated signal points, the lowest and largest values of the acceptance times efficiency (A × ϵ) for all three models are presented as a range. For the Z 0 B model, signals with DM mass m χ between 1 and 1000 GeV and mediator mass m Z 0 B between 1 and 2000 GeV are taken into consideration. The samples with the mediator mass m Z 0 ¼ 400-1400 GeV and pseudoscalar boson mass m A 0 ¼ 200-450 GeV are added for the Z 0 -2HDM model. For the heavy-scalar model, the values are taken from the signals points with m H ¼ 260 to 350 GeV and m χ ¼ 60 GeV. allowed. Pseudoscalar boson masses below 280 GeV are excluded for m Z 0 ¼ 1000 GeV and m χ ¼ 100 GeV. Tables VI and VII present the 95% C.L. observed and median expected limits on σðpp → hχχÞ × Bðh → γγÞ, respectively, for the Z 0 B benchmark model for different Z 0 B masses and the Z 0 -2HDM model for different pseudoscalar A 0 masses. The corresponding expected limits with one standard deviation are also shown.
A two-dimensional exclusion region in the plane formed by the mediator masses and the DM particle mass is obtained by means of a reweighting technique based on generator-level variables. Samples for a variety of mass points are generated using the MADGRAPH generator matched to PYTHIA 8.186 using the A14 tune for parton showering and hadronization. Bin-by-bin weights for each of the different mass points are obtained by comparing the generator level distribution of a given variable to the distribution for the same variable in a fully simulated benchmark sample. This procedure can be repeated for other variables. In the case of Z 0 B samples, weights are derived using the true E miss T . In the case of Z 0 -2HDM samples, weights are derived in a two-dimensional plane of the true E miss T and p γγ T . To validate this technique, several fully simulated and reconstructed samples are produced and their reconstruction-level variables are compared with the sample obtained from the reweighting technique. The acceptances of the samples after kinematic cuts agree within 5%, and the residual difference is treated as an extra systematic uncertainty in the signal yield. The observed and expected 95% C.L. limit contours for the signal strength σ obs =σ th are shown in Fig. 5 for both the Z 0 B and Z 0 -2HDM models, in which σ obs is the observed limit on the model cross section at a given point of the parameter space and σ th is the predicted cross section in the model at the same point. FIG. 4. Expected (dashed lines) and observed (solid lines) 95% C.L. upper limits on σðpp → hχχÞ × Bðh → γγÞ for (a) the Z 0 B model for g q ¼ 1=3, g χ ¼ 1, sin θ ¼ 0.3 and Dirac fermion DM m χ ¼ 1 GeV, and (b) the Z 0 -2HDM model for tan β ¼ 1, g Z 0 ¼ 0.8, m Z 0 ¼ 1000 GeV and Dirac fermion DM m χ ¼ 100 GeV, as a function of m Z 0 B and m A 0 , respectively. The masses of the neutral CP-even scalar (H 0 ) and the charged scalars (H AE ) from Z 0 -2HDM model are set to 300 GeV. The theoretical predictions of σðpp → hχχÞ × Bðh → γγÞ for these two models (dark-blue lines with blue bands representing their associated theoretical systematic uncertainties) are also shown. The inset shows a closeup view of the same figure in narrower ranges of both the x and y axes.  Figure 6 shows a comparison of the inferred limits to the constraints from direct detection experiments on the spinindependent (SI) DM-nucleon cross section in the context of the Z 0 B simplified model with vector couplings using the relation [5] in which μ Nχ ¼ m χ m N =ðm χ þ m N Þ is the reduced mass of the DM-nucleon system, and f p;n ¼ 3g q g χ =m 2 Z 0 B are the couplings between DM particles and protons and neutrons. The parameter Z is the number of protons in the considered nucleus and A the number of nucleons (both set to 1). Limits are shown at 90% C.L. For comparison, results from direct detection experiments (LUX [67], PandaX-II [68], XENON [69], superCDMS [70], and CRESST-II [71]) are also shown. The comparison is model dependent and solely valid in the context of this model. The results for the Z 0 B model, with couplings g q ¼ 1=3 and g χ ¼ 1 for this search, are more stringent than direct detection experiments for m χ < 2.5 GeV and extend to DM masses well below 1 GeV. The shape of the exclusion line at DM mass m χ ∼ 200 GeV to low masses is due to the loss of sensitivity in Z 0 B models where DM particles are produced via an offshell process. The impact of renormalization-group evolution effects [72,73] when comparing collider and direct detection limits is not taken into consideration here.  The results from this analysis, in which the region inside the contour is excluded, are compared with limits from the LUX [67], PandaX-II [68], XENON [69], superCDMS [70], and CRESST-II [71] 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 g q ¼ 1=3 and g χ ¼ 1. The impact of renormalization-group evolution effects [72,73] when comparing collider and direct detection limits is not taken into consideration here. C. Interpretation of the heavy-scalar model Figure 7 shows the m γγ distributions in the five categories as well as the fitted contribution of a heavy-scalar boson for illustration. No significant excess is observed in any category. In the heavy-scalar interpretation, the 95% C.L. upper limits on the σðpp→HÞ×BðH →γγχχÞ as a function of m H for m χ ¼ 60 GeV are shown in Fig. 8 and Table VIII, where a 100% branching fraction is assumed for H → hχχ. The upper limit at 95% C.L. is 15.4 fb for m H ¼ 260 GeV, and 4.3 fb for m H ¼ 350 GeV.

IX. SUMMARY
A search for dark matter in association with a Higgs boson decaying to two photons is presented. This study is based on data collected with the ATLAS detector, corresponding to an integrated luminosity of 36.1 fb −1 of proton-proton collisions at the LHC at a center-of-mass energy of 13 TeV in 2015 and 2016. No significant excess over the expected background is observed. For the Mono-Higgs category, a visible cross section larger than 0.19 fb is excluded at 95% C.L. for BSM physics processes producing missing transverse momentum and a SM Higgs boson decaying into two photons. Upper limits at 95% C.L. are also set on the production cross section times branching fraction of the Higgs boson decaying into two photons in association with missing transverse momentum in three different benchmark models: a Z 0 B model, a Z 0 -2HDM model and a heavy scalar boson (H) model. Limits at 95% C.L. are also set on the observed signal strength in a two-dimensional m χ -m Z 0 B plane for the Z 0 B model, and the m A 0 -m Z 0 plane for the Z 0 -2HDM model. Additionally, the results for the Z 0 B model are interpreted in terms of 90% C.L. limits on the dark-matter-nucleon scattering cross section, as a function of the dark-matterparticle mass, for a spin-independent scenario. For a dark-matter mass lower than 2.5 GeV, the constraint with couplings g q ¼ 1=3 and g χ ¼ 1 placed on the DMnucleon cross section is more stringent than limits from direct detection experiments. In the model involving the production of a heavy scalar boson, 95% C.L. upper limits are set on the production cross section times the branching fraction of H → hχχ → γγχχ for a dark-matter particle with mass of 60 GeV, where a 100% branching fraction is assumed for H → hχχ. The heavy-scalar model assuming H production through gluon-gluon fusion with a cross section identical to that of a SM Higgs boson of the same mass, is excluded for all the benchmark points investigated.

ACKNOWLEDGMENTS
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently. We acknowledge the support of ANPCyT FIG. 8. Observed and expected 95% C.L. upper limits on the σðpp → HÞ × BðH → γγχχÞ with a scalar DM candidate mass of 60 GeV as a function of the heavy-scalar-boson mass in the range 260 < m H < 350 GeV. A 100% branching fraction is assumed for H → hχχ. The theoretical prediction for the model (dark blue lines with blue bands representing their associated theoretical systematic uncertainties) is also shown. The theoretical cross section is assumed to be equal to that of a SM Higgs boson with the same mass produced in gluon-gluon fusion.