Observation of electroweak production of two jets in association with an isolated photon and missing transverse momentum, and search for a Higgs boson decaying into invisible particles at 13 TeV with the ATLAS detector

This paper presents a measurement of the electroweak production of two jets in association with a Z γ pair, with the Z boson decaying into two neutrinos. It also presents a search for invisible or partially invisible decays of a Higgs boson with a mass of 125 GeV produced through vector-boson fusion with a photon in the ﬁnal state. These results use data from LHC proton–proton collisions at √ s = 13 TeV collected with the ATLAS detector and corresponding to an integrated luminosity of 139 fb − 1 . The event signature, shared by all benchmark processes considered for the measurements and searches, is characterized by a signiﬁcant amount of unbalanced transverse momentum and a photon in the ﬁnal state, in addition to a pair of forward jets. Electroweak Z γ production in association with two jets is observed in this ﬁnal state with a signiﬁcance of 5.2 (5.1 expected) standard deviations. The measured ﬁducial cross-section for this process is 1 . 31 ± 0 . 29 fb. An observed (expected) upper limit of 0.37 (0 . 34 + 0 . 15 − 0 . 10 ) at 95% conﬁdence level is set on the branching ratio of a 125 GeV Higgs boson to invisible particles, assuming the Standard Model production cross-section. The signature is also interpreted in the context of decays of a Higgs boson into a photon and a dark photon. An observed (expected)


Introduction
Studying the self-couplings of the Standard Model (SM) vector bosons, precisely predicted through the SU(2) L ×U(1) Y gauge symmetry, provides a unique opportunity to better understand the electroweak sector of the SM and gain insight into possible anomalies due to new phenomena. Vectorboson scattering (VBS), V V → V V with V, V = W/Z / γ , is an interesting process to study, being sensitive to both the triple and quartic gauge-boson couplings, which are particularly sensitive to the presence of physics effects beyond the Standard Model (BSM) [1][2][3]. In hadron collisions, such as those produced at the Large Hadron Collider (LHC), VBS events occur whenever two vector bosons, radiated from the initial-state quarks, interact with each other [4], producing a final-state signature characterized by the two vector bosons and a pair of forward hadronic jets in opposite hemispheres. The two vector bosons can similarly annihilate and produce a SM Higgs boson in the so-called vector-boson fusion (VBF) mechanism. Depending on the subsequent decay of the SM Higgs boson, different final-state signatures can be exploited to investigate its properties. Precise knowledge of the various Higgs boson decay branching ratios is a fundamental aspect of understanding whether the 125 GeV scalar boson behaves according to the SM predictions or whether new physics phenomena modify the Higgs sector. This paper presents a set of SM measurements and searches for new phenomena in a final-state signature characterized by two forward hadronic jets, a photon, and a significant amount of unbalanced momentum in the plane transverse to the beam direction (E miss T ) due to undetected particles. In the SM this can result from V γ + jets production, where V is either a Z boson decaying into an undetected neutrino-antineutrino pair or a W boson decaying leptonically, where the charged lepton is not reconstructed in the detector. The latter is a background to the measurements and searches in this paper. The V γ + jets events are produced in the SM through a combination of 'strong' and 'electroweak' (EW) contributions: the former are produced through diagrams of order α 2 s α 3 at the Born level as shown in Fig. 1a, where α s is the strong coupling constant and α is the electromagnetic coupling constant, while the latter are produced more rarely through diagrams of order α 5 at the Born level 1 as shown in Fig. 1b, c. The Z γ + jets cross-sections have been computed at next-to-leading order (NLO) in α s for both the strong [5] and EW [6] production modes. It is not possible 1 One order of α is included for the decay of the vector boson.
to study VBS diagrams, as shown in Fig. 1b, independently of other electroweak processes (e.g. triboson production as shown in Fig. 1c) as only the ensemble is gauge invariant [7]. There is also interference between the SM electroweak and strong processes, which is accounted for in the measurement.
In Run 1 the Z γ + jets EW production cross-section has been measured in the dielectron and dimuon final states of the Z boson by the ATLAS and CMS experiments [8,9] with observed significances of 4.1 and 3.0 standard deviations, respectively. The CMS Collaboration measured EW Z (→ )γ jj after observing the process with a significance of 9.4 standard deviations in 137 fb −1 of 13 TeV proton-proton ( pp) collisions [10].
Profiting from the full 139 fb −1 dataset collected with the ATLAS detector during Run 2 of the LHC, the analysis presented in this paper reports the observation of EW Z (→ νν)γ jj production at the LHC.
The observation of EW production of Z (→ νν)γ jj lays the groundwork for further investigation of this signature in looking for possible hints of BSM physics, based on interesting and well-motivated benchmark scenarios involving new dark matter (DM) candidate particles or a hidden sector of new particles coupling with the SM Higgs boson. The existence of DM is evident from astrophysical observations [11], although its connection with the SM is still unknown because it has only been observed through gravitational interactions. In this paper, the connection of DM with SM particles is explored by introducing a coupling with the 125 GeV Higgs boson. A class of BSM scenarios, referred to as Higgs-portal models [12], feature a DM candidate behaving as a singlet under the SM gauge symmetries, with the Higgs boson playing the role of a mediator between the DM candidate and SM particles.
Two main benchmark models are chosen when searching for new phenomena in the considered final-state signature, probing the invisible decay of the Higgs boson or its semivisible decay into one invisible particle and one photon. The decay of the Higgs boson into invisible particles (H → inv.) when produced through the SM-predicted VBF process in association with an emitted photon [13,14] is probed in this paper for the first time. The Feynman diagram for this signature is shown in Fig. 2a. Compared to a similar signature without the photon, this one benefits from better background rejection and a higher signal reconstruction efficiency; however, it has a lower production cross-section.
No direct constraints on the H → inv. decays in this experimental signature currently exist. Constraints on the invisible Higgs boson branching ratio (B inv ) were set by the ATLAS Collaboration using the full Run 1 dataset at 7 and 8 TeV [15][16][17] and the full Run 2 dataset at 13 TeV [18,19]. Similar searches were carried out by the CMS Collaboration [20-22] using both the Run 1 and Run 2 datasets. The most stringent limits are from the statistical combination of the search results, for which ATLAS reports an observed (expected) upper limit on B inv of 0.26 (0.17) and CMS reports an upper limit of 0.19 (0.15), all at 95% confidence level (CL). In these combinations, the VBF production channel is the single channel with the highest expected sensitivity, for which ATLAS and CMS reported observed (expected) 95% CL limits on B inv of 0.37 (0.28) and 0.33 (0.25), respectively, using 36 fb −1 of Run 2 data.
In addition to the H → inv. benchmark, this analysis probes a dark-photon model [23-25] which predicts a light or massless 2 dark photon (γ d ) coupled with the Higgs boson through a U(1) unbroken dark sector. The Feynman diagram for this signature is shown in Fig. 2b, and BSM extensions to other resonance masses are possible [26]. In this model, a Higgs boson decays into a photon and an invisible dark photon with a branching ratio B(H → γ γ d ). CMS has probed this benchmark model by considering associated Z H production and VBF Higgs boson production [27,28], and reported observed (expected) 95% CL upper limits of 0.046 (0.036) and 0.035 (0.028) on B(H → γ γ d ) in their respective signa-tures. These two results were combined to yield an observed (expected) 95% CL upper limit of 0.029 (0.021).
The SM processes resulting in the same final-state signature as described before, mainly Z (→ νν)γ + jets and W (→ ν)γ + jets, represent background contributions in the searches for new phenomena. These processes were simulated, for both the Z γ + jets EW cross-section measurement and the searches for new phenomena, through dedicated Monte Carlo samples, detailed in Sect. 4. Minor and reducible background contributions arise due to the misreconstruction of objects in the detector, including hadronic jets reconstructed as photons. These background processes are estimated using simulation and data-driven techniques and are validated using dedicated data samples. The SM H → Z (→ νν)γ process produces the same signature as the investigated H → inv. + γ signal, but the small contribution of such events satisfying the criteria defined in Sect. 6 is neglected in the searches for new phenomena.
A brief description of the ATLAS detector is provided in Sect. 2. The dataset, simulated samples, and physics object selection are covered in Sects. 3, 4, and 5 , respectively. The analysis strategies and the event selection for the EW Z (→ νν)γ jj measurement and searches based on different signal hypotheses, H → inv. and H → γ γ d , are discussed in Sect. 6. The modelling of the different processes contributing to the considered signature is presented in Sect. 7, and the fit models and extracted results for the different signal hypotheses are described in Sect. 8.

ATLAS detector
The ATLAS experiment [29][30][31] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and near 4π coverage in solid angle. 3 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer. The inner tracking detector (ID) covers the pseudorapidity range |η| < 2.5. It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|η| < 1.7). The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to |η| = 4.9. The muon spectrometer (MS) surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering. A two-level trigger system is used to select events. The first-level (L1) trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz. This is followed by a software-based high-level trigger (HLT) [32,33] which reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions. An extensive software suite [34] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data samples
This analysis relies on the data collected by the ATLAS experiment from LHC pp collisions at √ s = 13 TeV during 2015-2018 stable beam conditions when all subdetectors were operational [35], corresponding to a total integrated luminosity of 139 fb −1 [36,37]. The data used in this analysis were recorded mainly using trigger algorithms based on the presence of missing transverse momentum, E miss T (described in Sect. 5) [38]. The trigger thresholds for the E miss T were determined by the data-taking conditions during the different periods, especially by the average number of multiple pp interactions in the same or neighbouring bunch crossings, referred to as pile-up. The first-level trigger threshold was 50-55 GeV, depending on the data-taking period, and the lowest-transverse-momentum ( p T ) unprescaled HLT thresh-Footnote3 continued axis along the beam pipe. The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of old for the E miss T trigger algorithm was 70 GeV in 2015, 90 GeV in 2016 and 110 GeV in 2017-2018. This corresponds to a trigger operating on its maximum-efficiency plateau for events with offline E miss T of ∼180 GeV, depending on the trigger thresholds, in this final state. Independent data samples exploited to study W (→ ν)γ + jets background were collected using the lowestp T unprescaled single-lepton triggers, with p T thresholds ranging from 20 to 26 GeV for the triggers with the tightest lepton identification criteria [39,40]. The E miss T triggers complemented the muon triggers to increase the number of accepted single-muon events by 28%, as discussed in Sect. 6.5.

Simulated event samples
Monte Carlo (MC) simulated samples are used to model both the SM and BSM processes. The full set of simulated samples is summarized in Table 1. The generated events were processed through a simulation [41] of the ATLAS detector geometry and response using Geant4 [42], and through the same reconstruction software as the collected data. For BSM signal samples with a Higgs boson mass different from 125 GeV, the detector response was simulated using a fast parameterized simulation of the ATLAS calorimeters [43] and the full Geant4 simulation for the other subdetectors.
The pile-up effects in the same and neighbouring protonbunch crossings were modelled by adding detector signals from simulated inelastic pp events to the original hard-scattering (HS) event. These were generated with Pythia 8.186 [44] using the NNPDF2.3lo set of parton distribution functions (PDFs) [45] and the A3 set of tuned parameters (tune) [46]. The energy scale and resolution for leptons and jets, their reconstruction and identification efficiencies, and the trigger efficiencies in the simulation are corrected to match those measured in data.

V γ + jets processes
The W (→ ν)γ + jets and Z γ + jets (together labelled V γ + jets) processes contributing to the signature considered in this analysis contain a charged lepton ( = e, μ or τ ) and a neutrino, a pair of neutrinos (νν) or a pair of charged leptons ( ) together with a photon and associated jets. The V γ + jets processes are split into two components based on the order in the electroweak coupling constant α. At tree level, the strong component is of order α 2 s α 3 and the EW component is of order α 5 ; example Feynman diagrams are shown in Fig. 1 for these different contributions. The strong component of each V γ + jets contribution was simulated for photon p T greater than 7 GeV using the Sherpa 2.2.8 [47] event generator at NLO precision in α s for up to one additional parton and at LO precision in α s for up to three additional partons. These cal-  [48] and OpenLoops [49] matrix element generators, and the parton-shower matching [50] was performed using the MEPS@LO [51] or MEPS@NLO [51][52][53][54] prescription. The NNPDF3.0nnlo set of PDFs [55] was used, along with dedicated parton shower tuning developed by the Sherpa authors. Electroweak radiative corrections to strong V γ + jets production have been computed at NLO [56][57][58] as weights in Sherpa, and these are roughly −2% to −4% relative corrections in the chosen signal region. A second sample of Z (→ νν)γ + jets was generated using Mad-Graph5_aMC@NLO [59] with the FxFx merging scheme [60] at NLO precision in α s for up to one additional parton, filtered for photon p T >10 GeV, and showered using Pythia 8 [61]. The difference between the Sherpa and MadGraph5_aMC@NLO Z (→ νν)γ + jets predictions is symmetrized around the Sherpa prediction and is taken as a source of modelling systematic uncertainty. Matrix elements for the EW contribution were calculated at LO in α s using MadGraph5_aMC@NLO 2.6.5, and the photon p T was required to be larger than 10 GeV. Generated events were showered using Pythia 8 [61] with the dipole recoil option enabled along with the CKKW-L [62,63] merging scheme. These EW samples are normalized to NLO QCD predictions obtained from VBFNLO [64] through a correction which depends on the dijet invariant mass (m jj ). The effect of the QCD factorization and renormalization scale choices is evaluated from the same NLO QCD process calculation in VBFNLO. The MadGraph5_aMC@NLO 2.6.5 EW samples include VBS contributions (shown in Fig. 1b) and triboson contributions (shown in Fig. 2c) which can depend on trilinear and quartic gauge couplings. The triboson processes contribute as much as 15% of the EW samples for low m jj values, 250 < m jj < 500 GeV, but less than 3% in the more signal-like high-m jj region. The interference between the EW and strong-production diagrams, which is of order α s α 4 , was simulated at LO in α s in Mad-Graph5_aMC@NLO 2.6.5 and the corresponding events were showered using Pythia 8 with the dipole recoil shower enabled. This simulated sample is not included in the background predictions but is used to calculate an uncertainty in the EW V γ + jets contribution, which is 5% or less in the m jj > 500 GeV kinematic region. Similarly, the EW production of a W or Z boson in association with two photons was found to be less than 2% of the EW V γ + jets samples calculated with MadGraph5_aMC@NLO 2.6.5, so it is not included.
The whole V γ + jets simulation (i.e. both EW and strong production) applies smooth-cone isolation [65] with a cone size of 0.1 and parameters n = 2 and = 0.10 to remove the collinear singularity between photons and charged partons which would otherwise appear in the process amplitude calculations.

V + jets processes
Simulation is used to model V + jets processes when the vector boson decays into neutrinos, muons or τ -leptons, and most of these simulated events are removed by the lepton vetoes, object overlap removal described in Sect. 5, and the removal of overlaps with V γ + jets, which is described later in this section. Therefore, their contribution is very small. However, V decays into electrons enter the selection mostly through the electrons being misidentified as photons, and this is modelled using the fully data-driven method described in Sect. 7.1. The strong-production V + jets sample was simulated with the Sherpa 2.2.7 event generator with the NNPDF3.0nnlo set of PDFs [55]. Parton-shower matching [50] was performed using either the MEPS@LO or MEPS@NLO prescription with associated parameters tuned by the Sherpa authors. The strong production of V + jets uses NLO matrix elements for up to two partons and LO matrix elements for up to four partons calculated with the Comix and OpenLoops libraries and the MEPS@NLO prescription. The samples are normalized to a next-to-next-toleading-order (NNLO) prediction for V + jets [66]. The EW V + jets sample was generated at NLO in α s , using Herwig 7 [67] to perform the parton shower and hadronization with the MMHT2014 2014nlo68cl PDF set [68]. Herwig 7 uses its Matchbox module [69] to assemble calculations with dipole shower algorithms [70] and interfaces with matrix element plug-ins from VBFNLO to compute NLO α s corrections in the VBF approximation.
After generating events for both the strong and EW V γ + jets processes as in Sect. 4.1, the samples may overlap with V + jets events in which an ISR/FSR photon is radiated. To avoid this overlap, V γ + jets events are removed if the photon passes the smooth-cone isolation and lies within R = 0.1 of an electron, muon, or τ -lepton that has p T > 10 GeV. In V + jets samples, events are accepted only if they satisfy the same criteria.

Top-quark processes
Minor background contributions originating from top-quarkrelated processes and associated production of top quarks with a W boson, were all modelled using the Powheg Box v2 [71][72][73] generator at NLO with the NNPDF3.0nlo PDF set. The events were interfaced with Pythia 8.230 for the parton shower and hadronization modelling with the A14 tune [74] and the NNPDF3.0nlo set of PDFs. The production of ttγ was modelled using the MadGraph5_aMC@NLO 2.2.3 generator at NLO. The events were interfaced with Pythia 8.186 for the parton shower and hadronization modelling with the A14 tune and the NNPDF2.3lo set of PDFs. The decays of bottom and charm hadrons were performed by EvtGen 1.6.0 [75] in all top-quark processes.

Additional background samples
The γ + jet background is simulated using Sherpa 2.2.2 at NLO in α s , similarly to the Sherpa V + jets samples in Sect. 4.1 except for the use of a dynamical merging scale to capture the fragmentation contribution [76]. Due to the large cross-section for this process, a partially data-driven 'jetsmearing' technique was used to increase the sample size as discussed in Sect. 7.1.
Other processes listed in Table 1 but not described in the previous subsections result in a negligible contribution to the analysis signature. Therefore, no further description is given.

Higgs boson processes
In the searches for H → inv. carried out in this paper the target is VBF Higgs boson production, although the small contribution from gluon-gluon fusion (ggF) Higgs boson production processes satisfying the selection criteria is also considered as a contributing signal. For both of these processes, corresponding MC samples were generated according to the details described in the following. The invisible Higgs boson decay was simulated using the SM H → Z Z * → 4ν decay with a 100% branching ratio. The difference in the relevant kinematic distributions between the H → 4ν process and Higgs boson decays into new undetected particles is negligible.
The VBF Higgs boson production process with an additional photon, as shown in Fig. 2a, was simulated for a Higgs boson mass of 125 GeV, requiring the presence of only electroweak vertices at LO and photon p T larger than 10 GeV. This process was computed to NLO accuracy in α s using MadGraph5_aMC@NLO interfaced with Herwig 7 [67,77] for parton shower and non-perturbative hadronization effects, using the PDF4LHC15 PDF set [78]. Parton shower uncertainties are computed from the relative difference between the process generated by Mad-Graph5_aMC@NLO at LO in α s with showering by Herwig 7 and the same process with showering by Pythia 8 [79] with the dipole recoil shower variable [80] turned on. The parton shower comparison is not made at NLO because the dipole recoil option in Pythia 8 is available only at LO, and without it there is a larger prediction of central emissions. 4 The same smooth-cone isolation as described in Sect. 4.1 is used to remove the collinear singularity between photons and charged partons. 4 The option gives an alternative approach to local recoils, where only one final-state parton takes the recoil of an emission. The dipole recoil option has been shown to better model the CMS data in VBF Z → processes [81,82]. When it is used, it improves the modelling of radiative emissions from non-colour-connected partons (the VBF jets), which are poorly modelled by Pythia 8 when it is turned off.
The ggF Higgs boson production process was simulated at NNLO accuracy in α s using Powheg NNLOPS [71][72][73]83,84], which achieves NNLO accuracy for arbitrary inclusive gg → H observables by reweighting the Higgs boson rapidity spectrum in MJ-MiNLO [85][86][87] to that in HNNLO [88]. The PDF4LHC15 PDF set [78] and the AZNLO tune of Pythia 8 were used. This simulation was interfaced with Pythia 8 for parton shower and nonperturbative hadronization effects. The ggF prediction from the MC samples is normalized to the next-to-NNLO crosssection in QCD plus electroweak corrections at NLO [89][90][91][92][93][94][95][96][97][98][99]. Photons present in these samples were generated using Photos [100,101] since no complete matrix element computation of ggF Higgs boson with a photon is available. 5 The simulation of ggF production of a 125 GeV Higgs boson described above is also used to model its decay into a photon and an invisible massless dark photon (γ d ) again normalized to a 100% branching ratio. The VBF production process was simulated to NLO precision in α s with the Powheg generator interfaced with Pythia 8 for hadronization and showering and the same γ γ d decay, as shown in Fig. 2b. The NLO electroweak corrections for VBF Higgs boson production were computed using HAWK [102] and were applied as a function of the Higgs boson's p T . The VBF samples were generated not only for a 125 GeV Higgs boson, but also for lighter and heavier Higgs bosons, for the interpretation of the search in the context of other scalar mediators. Throughout this paper, the samples with a Higgs boson assume SM couplings and use the narrow width approximation [89] for the various Higgs boson masses, ranging from 60 GeV to 2 TeV.

Object reconstruction
Objects are reconstructed from detector signatures by using identification algorithms widely deployed in ATLAS analyses. Candidate events are required to have a reconstructed vertex with at least two associated tracks, each with p T > 0.5 GeV and originating from the beam collision region in the x-y plane. The primary vertex in the event is selected as the vertex with the highest scalar sum of the squared p T of associated tracks [103].
Electrons are reconstructed by matching clustered energy deposits in the EM calorimeters to tracks in the ID [104], including the transition regions between the barrel and endcap EM calorimeters at 1.37 < |η| < 1.52. Electron candidates must have p T > 4.5 GeV and |η| < 2.47, and fulfill loose identification criteria. Depending on the p T and |η| range, muons are reconstructed by matching ID tracks to MS tracks or track segments, by matching ID tracks to a calorimeter energy deposit compatible with a minimum-ionizing particle, or by identifying MS tracks passing a loose requirement and compatible with originating from the IP [105]. Muon candidates must have p T > 4 GeV and |η| < 2.7, and fulfill a very loose identification criterion. No isolation requirement is placed on electron or muon candidates used as a lepton veto. In events with one or more leptons associated with the primary vertex, i.e. in control regions described in Sect. 6, muons (electrons) are required to satisfy medium (tight) identification criteria in order to improve the purity of background processes.
Photon candidates are reconstructed from clustered energy deposits in the EM calorimeter [104]. They must have p T > 15 GeV, fulfill tight identification and isolation criteria [104], and lie within |η| < 2.37 but not in the transition region (1.37 < |η| < 1.52) between barrel and endcap EM calorimeters.
Particle flow (PFlow) jets are reconstructed using the antik t algorithm [106,107] with a radius parameter of R = 0.4, using charged constituents associated with the primary vertex and neutral PFlow constituents as inputs [108]. The jet energy is calibrated with the effect of pile-up removed [109]. Jets are required to have p T > 20 GeV and |η| < 4.5. For jets with p T < 60 GeV and |η| < 2.5 the jet vertex tagger (JVT) discriminant [110] is used to identify jets originating from the HS interaction through the use of tracking and vertexing. The chosen JVT working point corresponds to a selection efficiency for HS jets of about 97%, evaluated on an inclusive Z (→ μμ) + jets sample. For |η| > 2.5, the two leading jets must pass a forward jet vertex tagger algorithm (fJVT) [111,112] that accepts 93% of jets from the HS interaction and rejects about 58% of pile-up jets with p T > 50 GeV, evaluated on an inclusive Z (→ μμ) + jets sample. Jets containing b-hadrons (b-jets) are identified using a multivariate discriminant (MV2c10) output distribution [113]. The working point is chosen to provide a 77% b-jet efficiency on an inclusive tt sample, with rejection factors of 6 and 134 for charm-hadron jets and light-flavour quark-or gluon-initiated jets, respectively.
To avoid double counting of energy deposits, the reconstructed objects are required to be separated according to the procedure detailed in Table 2. The overlap of photons with electrons, muons and jets is resolved by a R criterion. For leptons in the vicinity of jets, the R threshold for sufficient separation depends on the p T of the lepton to account for the collimation of boosted objects.
The unbalanced momentum in the transverse plane, referred to as missing transverse momentum or E miss T , is defined as the negative vectorial sum of the transverse momenta of all selected electrons, muons, photons, and jets, as well as tracks compatible with the primary vertex but not matched to any of those objects, this last contribution is the magnitude of the negative vectorial sum of all jets in the event with p T > 20 GeV before the JVT requirement. This is a powerful variable for rejecting events where the E miss T is generated by inefficiencies of the JVT requirement for HS jets.
Several cleaning requirements are applied to suppress non-collision backgrounds [116]. Misreconstructed jets can be caused by electronic noise, and jets from collisions are identified by requiring a good fit to the expected pulse shape for each constituent calorimeter cell. Beam-halo interactions with the LHC collimators are another source of misreconstructed jets. Those jets are identified by requirements on their energy distribution in the calorimeter and the fraction of their constituent tracks that originate from the primary vertex. The event is rejected if any selected jet is identified as a misreconstructed jet. Residual contributions of non-collision jets are absorbed into the normalization for the γ + jet background as described in Sect. 7.1.

Event selection
Based on the reconstructed objects in each event final state, the collected data are assigned to disjoint samples identified with specific phase-space regions, used in this analysis for different purposes. The signal regions (SRs) have the highest purity of signal process events, with mostly irreducible back-ground contributions. The control regions (CR) are enriched in background processes and are used to constrain estimates of such contributions. The validation regions (VR), similar in background process content to the CRs, are used to quantify the level of agreement between data and background predicted yields but are not included in the fit.
A baseline set of requirements is applied to select events with a photon, high E miss T , and the VBF jet signature considered in this paper, as described in Sect. 6.1. For the EW Z γ + jets cross-section measurement, events satisfying the baseline SR requirements are categorized according to their dijet invariant mass m jj , as described in Sect. 6.2. In the same section the fiducial volume considered for the crosssection measurement is defined. The SR-selected events for the H → inv. search are split into categories of different signal purities based on a multivariate analysis discriminant, which is described in Sect. 6.3. The SR selection and fitted discriminant are modified for the resonant kinematics of the H → γ γ d search, which are described in Sect. 6.4. In the statistical analysis of the data, described in Sect. 8, CR events are classified according to the same strategy described for the signal regions in the aforementioned measurement and searches.

Baseline event selection
Stringent background discrimination is made possible by the characteristic features of the VBF process, such as the presence of two highly energetic jets, typically in opposite η hemispheres of the detector, more forward than jets from non-VBF processes at comparable momentum transfer of the initial-state partons. These features lead, for example, to large values of the pseudorapidity separation η jj and invariant mass m jj of the jet pair. Multijet production predicted by QCD is instead characterized by two back-to-back leading jets in the transverse plane ( φ jj ∼ π ); therefore, to reduce that background contribution, the azimuthal separation of the leading jets, φ jj , is required to be smaller than 2.5 for all SRs.
The requirements defining the different regions considered in the analysis are summarized in Table 3. The following section describes the selection criteria.
The photon produced in association with either the Higgs boson (for the H → inv. search) or the Z boson (for the EW Z γ + jets cross-section measurement) is usually radiated from the scattering W bosons, and it is produced within the large rapidity gap between the two leading jets. For this reason the photon centrality C γ [117] is defined as where the subscripts 1 and 2 indicate the highest-and secondhighest p T jets in the event. The value of C γ is 1 when the photon is centred between the two jets characterizing the VBF signature,1/e when it is aligned with one of the two jets, and tends to zero when the photon is farther forward in |η| than the jets. Events are assigned to the SR if they satisfy a set of requirements that have been optimized to maximize the sensitivity of the analysis to the EW Z γ + jets and VBF H → inv. signals. These requirements are summarized in the following.
As discussed in Sect. 3, events are selected with the E miss Ttrigger algorithm. To ensure a trigger efficiency exceeding 97% in this topology, the offline E miss T , after full offline reconstruction and calibration of all the objects in the events, is required to be larger than 150 GeV. Corrections and uncertainties in those corrections to the simulation modelling of the E miss T trigger are discussed in Sect. 7.2. The leading and subleading jets are required to have p T > 60 GeV and 50 GeV, respectively, both satisfying the fJVT requirements mentioned in Sect. 5. These same two leading jets must be located in opposite η hemispheres (η( j 1 ) × η( j 2 ) < 0), to be well separated in pseudorapidity (| η jj | > 3.0) and not back-to-back in the plane transverse to the beamline ( φ jj < 2.5), and must have large invariant mass (m jj > 0.25 TeV).
Another characteristic of the EW processes in the VBF topology is reduced hadronic activity in the large rapidity gap between the two leading jets, caused by the absence of colour connection between the two quarks. To suppress the contribution from strong V γ + jets production with additional jets from QCD radiation in comparison with the benchmark signals, the equivalent centrality C 3 is defined for the third-leading jet in the event, if any, replacing η γ in Eq. (1) with the third-leading jet's pseudorapidity η 3 . A third jet with p T > 25 GeV can be present in the events, and additional jets are allowed only if they have p T < 25 GeV. The thirdleading jet is required to be in one of the two forward regions, corresponding to a small centrality value, C 3 < 0.7.
All jets must be well separated in azimuth from the E miss T direction, satisfying φ( j i , E miss T ) > 1 for each jet j i with p T > 25 GeV. At most one b-jet identified using the algorithm defined in Sect. 5 is allowed to be present, and the E jets,no-jvt T variable must be larger than 130 GeV in each event. The requirement on the maximum number of b-jets has a negligible effect (<0.1%) on the signal selection efficiencies, while it ensures that the dataset considered in the search for H → inv. produced through VBF presented in this paper is orthogonal to the dataset considered in the search for invisible decays of Higgs bosons produced in association with tt [118].
Each event must contain a single reconstructed photon with 15 < p T < 110 GeV, φ( E miss T , γ ) > 1.8, and C γ > 0.4. The upper bound on the photon p T reduces contributions from background processes, especially the γ + jet background. The z coordinate of the reconstructed photon extrapolated from the calorimeter to the beamline must be no more than 250 mm from the identified primary vertex to reduce the number of photons from non-collision backgrounds.
In the SR, all events with leptons are vetoed, specifically ensuring the absence of muon or electron candidates fulfilling the criteria described in Sect. 5. Electron and muon candidates considered in this veto are not required to satisfy any isolation or impact parameter requirement, and muons are considered before the overlap removal requirement.

Event classification and fiducial-volume definition for EW Z γ + jets measurement
In the signature considered in this analysis, the m jj observable helps in discriminating between the strong and EW Z γ + jets processes. The strong production Z γ + jets process tends to have a smaller m jj than the EW Z γ + jets process. In the context of the EW Z γ + jets cross-section measurement, the events satisfying the baseline SR requirements described in Sect. 6.1 are therefore binned according to their m jj value, with bin boundaries of 250 GeV, 500 GeV, 1000 GeV and 1500 GeV. The highest m jj in the selected data events is 4.0 TeV. The fiducial volume for the EW Z γ + jets cross-section measurement is defined to closely follow the SR definition reported in Sect. 6.1. Particle-level requirements are placed on stable particles (defined as having a mean lifetime cτ > 10 mm) before interaction with the detector but after hadronization. The four-momenta of prompt leptons, i.e. those which do not originate from the decay of a hadron, include the sum of the four-momenta of prompt photons within a cone of size R = 0.1 around the lepton, and the resulting leptons are referred to as 'dressed' leptons. Jets are constructed using an anti-k t algorithm with radius parameter of 0.4, excluding electrons, muons, neutrinos, and photons associated with the decay of a W or Z boson. Photon isolation momentum E cone20 T is defined as the transverse momentum of the system of stable particles within a cone of size R = 0.2 around the photon, excluding muons, neutrinos, and the photon itself. The photon isolation momentum is required to be less than 7% of the photon transverse momentum. The truth-E miss T is defined as the vector sum of the transverse momenta of all neutrinos in the simulated final state, including those resulting from hadronization. A veto on dressed muons or electrons with p T > 4 GeV and |η| < 2.47 is applied. Jets and photons are required to be separated in angular distance, with R( j, ) > 0.4. The remaining requirements, similar to the SR ones, are summarized in Table 4, and the acceptance times reconstruction efficiency for the fiducial-volume definition is 0.4%.

Event classification for H → inv. search
In the context of the search for H → inv. signal, a Dense Neural Network (DNN), based on the Keras library [119, 120] with Tensorflow [121] as a backend, was designed and trained to separate such signal from the background contributions by using kinematic features. The DNN assigns signal-like (background-like) events an output score close to 1 (0).
The DNN architecture is composed of three blocks, each of which includes a Dense layer with 384 neurons, a Dropout rate and a Batch Normalization. The Rectified Linear Unit (ReLU) activation function [122] is used for each Dense layer, while a sigmoid activation function is used for the output node. To regularize the training and avoid overtraining, the Dropout and Batch Normalization layers are included in the DNN architecture. Each of the three blocks receives as input the output of the previous block concatenated with the input features.
The DNN is trained using a mixed sample of signal and background events, weighted according to their expected yields in the considered final state. To increase the number of events available for the training, the selection criteria defining the training sample phase space are looser than those for the SR and CRs. Events used for the training are required to have either two or three jets satisfying the identification criteria described in Sect. 5, a single photon, and E miss,lep-rm T > 140 GeV. The sample for each considered process is split into a training sample and a testing sample composed of 85% and 15% of the original sample, respectively. The training sample is further divided into a subsample used in the proper training (80%) and a sub-sample for validation (20%).
During the DNN training the weights are updated via an optimizer used to minimize the binary Cross-Entropy loss, while at each step a customized metric is evaluated, representative of the expected limit E on the H → inv. signal. This expected-limit approximation is defined below and the expected-background error in the numerator includes the background yield, the MC statistical uncertainties, and the data background uncertainty in data as provided by the CR: is the contribution to the background from W γ + jets and Z γ + jets processes, d i (CR) is the number of events in data in the onelepton CR defined in Sect. 6.5, and w i, j is the MC weight of each event (ev) normalized to the integrated luminosity summed over the events in the considered bin.
After each training epoch the DNN weights are stored only if E is improved and the difference in loss between the training and validation samples is smaller than 3%, which prevents failures from overtraining.
After the composition of the training samples and the DNN architecture were established, a backward elimination procedure was carried out to reduce the number of input features with minimal performance loss. Input features are varied up and down by 10% of their nominal values while evaluating the impact of each as the change in the mean DNN output score. The most relevant input features are expected to have a greater impact on the output if their value is modified, while the least relevant ones will produce smaller variations, and hence can be dropped with no significant loss. When an input feature is dropped the DNN is retrained with the new subset of input features, and the procedure is repeated while the relative change in E from its value with that with all input features is less than 2%. The remaining features have a non-negligible impact on the DNN; hence they are used in the analysis. The 8 most significant kinematic features were selected: the separation of two leading jets in η, η jj , and in φ, φ jj ; the dijet invariant mass m jj ; the photon pseudorapidity η(γ ); the two leading jet transverse momenta p j 1 T and p j 2 T ; the subleading jet pseudorapidity η( j 2 ); and the leptonsubtracted missing transverse momentum E miss,lep-rm T . When their nominal value is perturbed by 10%, η jj , E miss,lep-rm T and m jj produce the largest change in DNN output score; therefore, they are interpreted as the most important in separating signal from background events.
All events are categorized according to DNN output score into four bins chosen to optimize E, subject to the requirements that each bin be at least 0.05 units wide in output score and that there be at least 20 background events expected in each bin in order to have sufficient background yield in the CRs to validate the predictions. The optimal boundaries of the DNN output score bins are [0.0, 0.25, 0.6, 0.8, 1.0]. The acceptance times reconstruction efficiency for VBF Higgs boson production with an additional photon is around 0.7%, which comes from the high E miss T threshold as well as the other selections.

Event classification for H → γ γ d search
The H → γ γ d decay presents a striking signature characterized by a kinematic edge at the Higgs boson mass in the distributions of the transverse mass constructed from the photon and E miss T , defined as: To increase the search's sensitivity to this semi-visible signal, a dedicated signal region SR γ d is defined by changing some of the selection requirements in Sect. 6.1. No requirement on φ( E miss T , γ ) is applied when targeting this specific benchmark and, for m T (γ , E miss T ) > 150 GeV, the upper bound on the photon p T is relaxed to 0.733 × m T (γ , E miss T ). The requirement on φ jj is tightened, so that all events have φ jj < 2.0. These selections are changed for both the SRs and the CRs defined later in Sect. 6.5. The events satisfying the SR γ d requirements are therefore binned according to m T (γ , E miss T ) value, with bin boundaries of 0 GeV, 90 GeV, 130 GeV, 200 GeV and 350 GeV. The events are further separated according to the value of the dijet invariant mass m jj , into events with m jj < 1.0 TeV and events with m jj ≥ 1.0 TeV. These two m jj regions have different relative contributions of the VBF-and ggF-produced Higgs boson signals, so the separation improves the overall sensitivity to the H → γ γ d decay. The acceptance times reconstruction efficiency for VBF H → γ γ d production is 0.4% for a 125 GeV scalar mediator but increases to 3.1% for a 500 GeV scalar mediator.

Control region definitions
Mutually exclusive CRs are defined in order to constrain the normalization and validate the modelling of kinematic distributions for the two dominant backgrounds from W (→ ν)γ + jets and strong Z (→ νν)γ + jets production. These CRs have kinematic features similar to those in the SR but are dominated by the background processes. The EW production of Z (→ νν)γ jj events, which is measured in this paper and constitutes an irreducible background contribution in the search for invisible and semi-visible decays of the Higgs boson, is mitigated through the event classification applied to SR events described in Sects. 6.3 and 6.4.
The W γ μν CR and W γ eν CR are defined by changing the zero-lepton requirement for the SR to instead require exactly one lepton, either a muon or an electron, respectively. Data events in the region characterized by one electron are collected using a single-electron trigger algorithm, while those characterized by one muon are collected using either the same E miss T trigger algorithm used to select SR events or a singlemuon trigger algorithm. In both cases the leptons are required to have p T > 30 GeV in order to be on the single-lepton trigger efficiency plateau.
In each of these CRs the presence of a single muon (electron) satisfying the medium (tight) identification criteria, passing the loose isolation requirement, and surviving the overlap removal procedure is required. The efficiency of selecting muons (electrons) averaged over p T > 30 GeV and |η| < 2.7 (2.47) is larger than 96% (80%). To ensure that electron and muon candidates originate from the primary vertex, the significance of the track's transverse impact parameter relative to the beam line must satisfy |d 0 /σ (d 0 )| < 3 (5) for muons (electrons), and the longitudinal impact parameter z 0 is required to satisfy |z 0 sin(θ )| < 0.5 mm. As done in the SR, the looser identification criteria are applied to define the veto on any additional lepton.
To reduce the contribution from jets faking electrons, E miss T > 80 GeV is required. The Fake-e CR, containing events with E miss T < 80 GeV, is used to estimate the jetsfaking-electrons background. The requirements on the reconstructed jets and photon are the same as for the SR, while requirements on the E miss T , other than the one specifically for the Fake-e CR mentioned above, are replaced by equivalent ones on the E miss,lep-rm T .
To check the modelling of the strong Z (→ νν)γ + jets background contribution, a Z γ Rev.Cen. CR is defined by applying all the requirements of the SR definition but reversing the photon centrality requirement, selecting events with C γ < 0.4. There is a 7% EW Z (→ νν)γ jj contribution to the Z γ Rev.Cen. CR, which is taken into account in the statistical analysis of the cross-section measurement for this process.

Validation region for γ + jet background
In order to correct the normalization of the γ + jet background process and validate the modelling of its kinematic features, an orthogonal sample of events characterized by a small amount of E miss T in the final state is considered. In this region (Low-E miss T VR) the same selections as the SR are applied, except that the E miss T has to be in the range 110 < E miss T < 150 GeV, the upper bound on the photon p T is removed, no E jets,no-jvt T requirement is applied, and the leading dijet system's invariant mass has to satisfy 0.25 < m jj < 1.0 TeV.
For the H → γ γ d search, a similar Low-E miss T validation region is defined (Low E miss T -γ d VR) with the same selection as the Low-E miss T VR but also requiring m T (γ , E miss T ) < 110 GeV without any selection on φ( E miss T , γ ).

Data analysis
The different analyses presented in this paper focus on processes sharing a similar signature. For the EW Z γ + jets cross-section measurement, all the SM processes contributing to the considered final-state signature are estimated using both the simulated samples, described in Sect. 4, and datadriven techniques, described in Sect. 7.1. The m jj distribution is used to categorize events into SR selections with different signal-to-background ratios, as detailed in Sect. 6.2. In the search for invisible Higgs boson decays, the same strategy is adopted to estimate the SM backgrounds but a different observable is considered. A multivariate discriminant is designed to categorize events falling in the SR according to an increasing expected significance for the invisible Higgs boson decay signal contribution relative to the SM background expectation, as detailed in Sect. 6.3. For the H → γ γ d search, events are instead classified according to their kinematic features, as detailed in Sect. 6.4. Control regions enriched in background contributions, described in Sect. 6.5, are considered in the statistical analysis of the data, while validation regions are considered only to qualitatively assess the modelling of specific background contributions, as detailed in Sect. 6.6. The theoretical and experimental uncertainties of the signal and background predictions are described in Sect. 7.2.

Background contribution estimation
The dominant contributions entering the analysis SR are Z (→ νν)γ + jets and W (→ ν)γ + jets events in which the lepton from the W decay is lost mostly because it falls outside of the p T or η acceptance. Roughly 28% of W γ + jets events in the 0-lepton SR come from hadronic τ -lepton decays with truth p T > 20 GeV and |η| < 2.5; this fraction is too small to benefit from a veto on reconstructed hadronic τ -lepton decays, which would incur additional uncertainties on their reconstruction efficiency. The uncertainty in the modelling of the lost lepton is included in the MC simulation, as detailed in Sect. 7.2. In order to validate the theoretical and experimental uncertainties, these major backgrounds are tested in the CRs described in Sect. 6 characterized by the same requirements as those of the SR but requiring events with exactly one lepton or reversing the C γ requirement. The CRs with one selected lepton (e.g. W γ eν CR and W γ μν CR) are expected to be dominated by W (→ ν)γ + jets events while the Z γ Rev.Cen. CR with the reversed C γ requirement is expected to be dominated by strong production of Z (→ νν)γ + jets events. The background yield compared with the data in the aforementioned CRs is utilized to determine the normalization of the W (→ ν)γ + jets and Z (→ νν)γ + jets backgrounds by allowing their normalization to float in the different fit models, as described in Sect. 8.
An additional background contribution results from Z + jets or W + jets events in which one of the jets is misreconstructed as a photon. The simulation described in Sect. 4 is not expected to properly reproduce the rate of jets misreconstructed as photons. Therefore, the jet → γ misreconstruction rate is estimated using a technique tuned on data and similar to the one described in Ref. [123]. The method, often referred to as an 'ABCD' method, relies on the evaluation of the photon candidate event yield in a two-dimensional plane defined by the amount of transverse energy deposited in clusters of cells within a cone of size R = 0.4 around the photon, excluding the photon cluster itself (i.e. isolation), and by the quality of the photon identification criteria (tightness) [104]. A photon signal region (region A) is defined by photon candidates that are isolated and satisfy the tight identification requirements described in Sect. 5. Three background regions are defined in the isolation-tightness plane, consisting of photon candidates which are tight but non-isolated (region B), non-tight but isolated (region C) or neither tight nor isolated (region D). A photon candidate is defined as nonisolated if it does not satisfy the requirement on the amount of isolation transverse energy, while it is classified as nontight if it fails the tight identification but satisfies a modified set of requirements related to four of the selections associated with the shower-shape variables computed from the energy deposits in the first layer of the EM calorimeter. No correlation between the photon identification and isolation requirements would imply that the numbers of photon candidates in the four regions ( A, B, C, D) satisfy the condition N A /N B = N C /N D . Although this condition is almost fully satisfied, the residual correlation between the photon isolation and tightness is taken into account through correction factors estimated from MC simulations, differing from unity by up to 40%. A further correction to the described method is included to take into account the contamination from real photons in the three regions which are supposed to be dominated by fake photons (B, C, D). This contribution is evaluated using MC simulation and is parameterized through coefficients representing the fraction of real photons in each of the aforementioned regions relative to the fraction of real photons observed in region A. With these two corrections, this method determines the jet → γ misreconstruction rate and consequently the contribution due to this background in all the kinematic regions considered in this analysis. The nontight photon and isolation definitions were varied to compute a relative systematic uncertainty, estimated to be 80-90% for this background process.
The background contribution due to electrons misidentified as photons is estimated through an e ± → γ misreconstruction rate (fake rate), determined from a comparison of the rates of Z boson reconstruction in the e ± γ and e + e − final states, as in Ref. [124]. This background is small in the EW Z γ + jets measurement and in the Higgs boson invisible decays interpretation but is more relevant in the dark-photon interpretation. The full Run 2 dataset is used to select Z → ee events in which the electron pairs in the final state are either reconstructed as an e + e − pair or misreconstructed as an e ± γ pair. The invariant mass m ee or m eγ is required to be consistent with the Z boson mass within 10 GeV. The yields of Z events are then obtained from a simultaneous fit to the m ee or m eγ distributions, in order to subtract the contamination from jets misidentified as electrons or photons in the two samples using an extrapolation from the sidebands of the pair's mass. The e ± → γ misidentification rate, measured as a function of |η| and p T , varies between 1.5% and 9%, with larger values associated with larger values of |η|, since the misidentification rate depends on the amount of material in front of the calorimeter. The background contribution from misidentified electrons is evaluated in the different regions considered in this analysis by applying the calculated e ± → γ misidentification rates to event yields resulting from the same criteria except that an electron is required instead of a photon. The misidentification rate includes uncertainties from varying the fitted range of the pair mass, turning the jet background subtraction on and off, and accounting for biases in the energy of the photon compared to that of the electron. The total uncertainties vary from 30% for p T = 15 GeV to 15% for p T = 100 GeV.
The γ + jet background is a minor one in this analysis since these events do not have intrinsic E miss T , other than the small contribution from the presence of neutrinos in heavy-flavour jet decays. While only a small fraction of such events exceed the E miss T requirement of the SR definition, most of this background contribution results from jet mismeasurement, yielding a significant amount of misreconstructed E miss T . The sample of simulated events for this large cross-section process is limited, so a jet smearing approach is used to increase the sample size by a factor of 20. A simulated sample of γ + jet events is filtered with a set of criteria looser than those used in the SR definition so that events include one identified photon and two or more reconstructed jets with | η jj | > 2.5. The leading and subleading reconstructed jet p T must be larger than 50 and 30 GeV, respectively. The energy of each truth jet, which includes neutrino momenta, in these events is smeared according to the jet energy response evaluated in bins of p T of the corresponding truth-level jets; jet φ and η are smeared with a Gaussian distribution according to the angular resolution of PFlow jets as measured in 13 TeV data [125]. The pile-up jet energies are not smeared, but they are included in the jet counting. After all the truth jets have been smeared, the corresponding JVT and fJVT tagging scores are recalculated for each jet. The true-E miss T magnitude and direction due to heavy-flavour decay neutrinos are resampled from simulated sample truth information. In addition to this redefinition of the E miss T in the event, corrections due to jet smearing are also propagated to the reconstructed E miss T as well as all other kinematic quantities used to define the SR. The γ + jet sample obtained with this procedure is normalized according to the simulation-to-data comparison in the Low-E miss T VR described in Sect. 6, after subtracting all the other background contributions in that region. Normalization to data absorbs any residual contributions of detector noise or other misreconstructed jets. The resulting normalization, which was found to be statistically consistent over all data-taking years, is 0.91 ± 0.36 with the uncertainty coming from the statistical uncertainties of the data. Uncertainties in the extrapolation from the Low-E miss T VR are computed using the difference between the E miss T trigger efficiencies for γ + jet and Z (→ νν)γ + jets events, which results in a 75% uncertainty in yield. Lastly, the smearing parameterizations are varied within their uncertainties, and the truth-E miss T resampling is turned on and off. The total uncertainty in this background is roughly 90%.
Some of the events in the W (→ ν)γ + jets CR are events where a jet or a photon is misidentified as a lepton. Events of this kind result from γ + jet production in which a jet is misidentified as a lepton, diphoton production where a photon mimics the prompt electron, or to multijet events where one jet mimics the photon and another one mimics the lepton. The contribution from fake-lepton events is found to be negligible in the case of events with one muon. To evaluate this contribution to the W γ eν CR, a corresponding CR (Fake-e CR) is defined, as discussed in Sect. 6.5. The ratio of the yields of fake electrons in the W γ eν CR and Fake-e CR, labelled R fake-e , is estimated in the same regions except that the one electron satisfying the tight identification requirements is replaced by an electron satisfying a loose identification requirement but failing the tight W γ eν CR electron definition. This is called the anti-ID selection. The number of data events in the anti-ID W γ eν CR, after subtracting residual contributions from prompt-lepton events using simulation, is divided by the number of data events in the anti-ID Fake-e CR to compute R fake-e . The ratio R fake-e is measured to be 0.14 ± 0.11 in the baseline selection (Sect. 6.1) and 0.26 ± 0.10 in the selection used in the H → γ γ d search (Sect. 6.4), with the uncertainties coming from statistical uncertainties of the data. This ratio is used to scale the fake-electron contributions in the Fake-e CR to obtain the expected background in the W γ eν CR.

Systematic uncertainties
Theoretical uncertainties affect all simulated signal and background processes. They originate from the limited order at which the matrix elements are calculated, the matching of those calculations to parton showers, and the uncertainty of the proton PDFs.
For the Z γ + jets and W γ + jets processes the higherorder matrix element effects and parton-shower-matching uncertainties are assessed by varying the renormalization and factorization scale choices used in the event generation. The effect of resummation and CKKW matching scale [63] choices are found to be negligible. The factorization and renormalization scales are varied up and down by a factor of two using event weights varied in the Sherpa MC simulation of strong production, while they are varied in separate calculations in VBFNLO for EW production. The corresponding uncertainties are calculated by taking an envelope of the seven factorization/renormalization scale variations: the central value, each scale independently varied up/down, and both scales coherently varied up/down. For the strong-production V γ + jets background, the effect of the 7-point scale variations on the expected yield ranges from ∼25% to ∼56% in the kinematic regions considered; the corresponding values for the EW V γ + jets process range from ∼3% to ∼11%. Since the EW V γ + jets simulation is leading-order Mad-Graph5_aMC@NLO reweighted by NLOVBFNLO, the 7point scale variations are computed from VBFNLO. For the strong-production V γ + jets background the scale variations are propagated to both the matrix element and the parton shower, while for the EW V γ + jets samples the Pythia 8 A14 'eigentune' variations [74] are summed in quadrature to assess the uncertainty in the parton shower modelling. The systematic uncertainty is found to be in the range 4-15%.
For the Z (→ νν)γ + jets process, the interference between EW and strong production is computed at LO in MadGraph5_aMC@NLO + Pythia 8, as discussed in Sect. 4.1. The full prediction of the interference sample is assigned to the EW Z (→ νν)γ jj sample as a one-sided uncertainty, as is done in Ref. [8], in all SR and CR bins, and its largest contribution is at large m T (γ , E miss T ) with a −22% uncertainty in the prediction without interference.
An uncertainty due to the choice of MC generator is evaluated by comparing strong V γ + jets samples produced with Sherpa and MadGraph5_aMC@NLO: the full difference between Sherpa and MadGraph5_aMC@NLO is considered as a systematic uncertainty, and the difference in the predicted yields is symmetrized around the Sherpapredicted central value for this background yield. The difference is as large as 20% depending on the signal region or control region bin, but the biggest differences are at large m T (γ , E miss T ). Both the shape and normalization differences are accounted for in the fit model.
The PDF uncertainties of the W γ + jets and Z γ + jets backgrounds are evaluated in each region and each bin of the SR as the standard deviation of the 100 PDF replicas of the NNPDF set. The overall uncertainty is found to be smaller than a few percent.
For the Higgs → γ γ d signal samples, the VBF and ggF Higgs boson inclusive production cross-sections and uncertainties are provided by the LHC Higgs Cross Section Working Group [89]. For the VBF process, m jj -dependent renormalization and factorization uncertainties and their correlations were computed by the LHC Higgs working group [89] for the Powheg+Pythia 8 signal samples. The parton shower uncertainty for the VBF signal is estimated to range from 2% to 4% by comparison with a Powheg+Herwig 7 sample. There is a 2% uncertainty due to the Higgs boson p Tdependent NLO electroweak correction, which is estimated from HAWK [102]. Uncertainties due to the PDF choice are assessed from the average variations of the PDF4LHC15 set of PDFs. For the ggF process, uncertainties from the renormalization and factorization scale variations are also assigned in categories of Higgs boson p T , number of jets, and m jj by the LHC Higgs Cross Section Working Group [89]. The resulting uncertainties are 20-30%. The smaller PDF and parton shower variations are also included. The same uncertainty treatment is used for the small contribution of ggFproduced H → inv. events.
For the VBF Higgs+γ sample, the factorization and renormalization scales are varied up and down by a factor of two through event weights varied in the MC matrix element and parton shower simultaneously. The envelope of the seven factorization/renormalization scale variations considered is assigned as the corresponding 1-2% systematic uncertainty. The uncertainty due to the modelling of the parton showering is assessed by comparing relevant kinematic distributions for Pythia 8 and Herwig 7 samples generated at LO in α s , as discussed in Sect. 4.5. The relative differences of 2-8% are applied to the NLO signal sample as uncertainties. Uncertainties of 1-2% due to the PDF choice are assessed from the average variations of the PDF4LHC15 set of PDFs.
Several experimental uncertainties impact the sensitivity of the analyses presented in this paper. They are grouped into categories: uncertainties in the luminosity, uncertainties in the trigger efficiencies, and uncertainties related to the reconstruction of physics objects such as electrons, muons, jets, and the E miss T . A summary of the systematic impact on the measured signal strength or limits set by this search is reported in Table 7.
The uncertainty in the luminosity is 1.7% [36] and impacts the signal yield and the simulated background yield.
The E miss T requirement of 150 GeV in the SR is key to maximizing the sensitivity of the analysis. The E miss T triggers are not fully efficient, so systematic uncertainties are used to account for possible trigger efficiency differences between data and simulation. This is done by comparing the combined L1+HLT trigger efficiency as a function of E T -trigger efficiencies for simulated W (→ μν) + jets, W (→ μν)γ + jets, and Z (→ νν)γ + jets are found to be statistically consistent in all the considered regions, so the same correction for differences between the data and simulation trigger efficiencies is applied to all simulated events passing the E miss T trigger. The correction as a function of E miss T varies from around (3-6)% ± 4% at E miss T = 150 GeV to less than 0.4% ± 1% for E miss T > 200 GeV depending on the specific E miss T -trigger algorithm used in a given data-taking period. The uncertainty in the correction comes from the data statistical uncertainties used to derive the correction. For the W γ ν CR scale factors and uncertainties in events passing lepton triggers, the corresponding single-lepton triggers corrections are applied [39,40].
Systematic uncertainties are calculated for lepton reconstruction and isolation efficiencies [104,105] and for the energy scale and resolution [104]. For the electron (muon) veto, an uncertainty in the electron (muon) reconstruction inefficiency is taken into account. For jets, uncertainties are derived in the energy scale and resolution [109] and for the pile-up tagging efficiencies [39,110]. For the photon, uncertainties in the reconstruction, isolation, energy scale and resolution [104,126] are considered. The above uncertainties associated with the reconstructed objects are propagated to the calculation of E miss T , as is the uncertainty in the scale and resolution of the E miss T soft term [115].

Fit models and results
The statistical analysis carried out has two objectives: first measuring the EW Z γ + jets production cross-section, and then searching for evidence of BSM physics in specific models involving invisible or semi-visible decays of the Higgs boson. The different processes of interest are compared to data using a profile-likelihood-ratio test statistic in a frequentist approach. A maximum-likelihood fit to the observed data in each bin of the event classification described in Sects. 6.2, 6.3, and 6.4 is used to set constraints on the signal strength μ for each model considered, all using asymptotic formulae [127]: a two-sided confidence level (CL) interval with the CL s definition [128] is extracted for the EW Z γ + jets cross-section measurement, while one-sided confidence levels calculated with the same approach are considered to set upper limits on new physics contributions.
Each bin i of the SR is assumed to include a number of events N i where the total expected background yield is B i and the signal contribution is S i . A likelihood function L is defined as where P(x|ν) is the Poisson probability density function, G(x|θ) is the probability density function of a Gaussian with unit width, and θ j represents the nuisance parameters corresponding to each considered uncertainty. The expected background yield in a given bin, B i , is given by the sum of several contributions, where a normalization factor β for each background component is included. The β factors represent the overall normalization of a given background contribution estimate. In the maximum-likelihood fit such normalization factors can be either fixed to 1 or left floating in the maximization. The treatment of these factors in the statistical analysis framework is detailed in the following for each scenario considered. In the framework of the EW Z γ + jets crosssection measurement, this process is considered as the signal and its corresponding β factor is replaced by the parameter of interest μ Z γ EW . In the framework of the searches for invisible or semi-visible decays of the Higgs boson, expected signal yields are evaluated with the assumption B inv = 1 or B(H → γ γ d ) = 1, which allows the parameter of interest μ to directly determine the actual Higgs boson decay branching fraction to the considered invisible or semi-visible particles. 6 The signal considered in this framework is normalized to the SM cross-section for Higgs boson production. The expected yields depend not only on B i and S i but also on the nuisance parameters, although these parameters are constrained in the likelihood function by the G(0|θ j ) factors. Some of the systematic uncertainties affecting background and signal predictions vary the final observable by less than 0.1% and these are ignored to improve the stability of the fit with no loss of accuracy. Each experimental uncertainty source is taken to be fully correlated across all signal and control regions. The applied correlation of theory systematic uncertainties depends on the type of uncertainty. PDF uncertainties are treated as fully correlated across bins. The uncertainty in the interference between the EW and strong Z (→ νν)γ + jets processes is treated as a fully correlated one-sided uncertainty in the EW Z (→ νν)γ jj process. The parton showering in each of the four V γ + jets samples (EW and strong-production components of W γ + jets and Z γ + jets) is correlated across all bins and similar for the scale variations in the strong-production V γ + jets samples, except the scale uncertainties in the onelepton CRs which are not correlated with the other bins. A separate nuisance parameter for the latter is motivated by the less than 1% background contribution to the one-lepton CR, and the contribution due to missing one of the leptons. Scale uncertainties for EW Z γ jj and EW W γ jj are treated as uncorrelated among regions and correlated in bins of the same region to avoid constraining this source of uncertainty. Within rounding, no change in the measured significance nor in the extracted limits was found compared with correlating these uncertainties across all analysis bins.

Fit model and results for the EW Z γ + jets cross-section measurement
The signature considered in this paper has a significant contribution from EW Z (→ νν)γ jj production, which has not been observed previously. A measurement of the cross-section for this process is an important first step, as a result complementary to the prior SM CMS observation of EW Z (→ )γ jj [10]. No BSM signal contributions are considered in this subsection.
For this measurement events are categorized into four m jj bins, as described in Sect. 6.2, to profit from differences between the kinematic distributions of EW Z γ + jets production and strong Z γ + jets production. The statistical analysis of data entering either the SR, W γ μν CR, W γ eν CR, Z γ Rev.Cen. CR, or Fake-e CR is performed according to their m jj bins, expanding the likelihood function defined in Eq. (2) as The EW Z γ + jets signal contribution entering the Z γ Rev.Cen. CR is taken into account in the statistical analysis. The fake-electron background in the CR bin at high E miss T is determined by the product of B fake-e and the transfer factor R fake-e (see Sect. 7.1). The symbol B Fake-e CR fake-e represents the number of events with a fake electron in the corresponding bin of the Fake-e CR at low E miss T . For this measurement, in addition to the parameter of interest μ Z γ EW , the normalization factors β Z γ strong and β W γ corresponding to the strong-production Z γ + jets component and, inclusively, to the W γ + jets component of the investigated distribution respectively, are allowed to float in the fit.
The result of the maximum-likelihood fit to data in the 4 SR bins and 16 CR bins is shown in Fig. 3, with the bestfit model propagated in all the regions. The SR bin-by-bin yields and CR yields are shown in Table 5 for the SM process contributions and data. The level of agreement between the data and the prediction from background simulation is good overall and is better after the fit, as shown in the lower panel of Fig. 3.
The best-fit values of the three free-floating normalization factors appearing in Eq. (3) are reported in Table 6. In particular, the EW Z γ + jets normalization best-fit value is 1.03 ± 0.16(stat) ± 0.19(syst) ± 0.02(lumi). The excess over the background-only hypothesis is quantified by a p-value using the profile likelihood ratio, evaluated at μ Z γ EW = 0, as a test statistic. EW Z γ jj is observed with a significance of 5.2σ with respect to the other SM background processes, and the expected significance is 5.1σ . The statistical component of the uncertainty includes only the data statistical uncertainty with other sources such as the limited number of   Table 6 The best-fit values and corresponding uncertainties of the three free-floating normalization factors derived from the statistical analysis described by Eq. (3) 1.03 ± 0.25 1.02 ± 0.41 1.01 ± 0.20 simulated events and the normalization of the backgrounds included in the systematic component of the uncertainty. The impact on the measurement of μ Z γ EW from different groups of uncertainties is shown in Table 7. It is evaluated by repeating the fit procedure, after fixing the nuisance parameters corresponding to each group of systematic uncertainties, in turn, to their best-fit values, and subtracting the new vari-ance (σ 2 ) of the best-fit value of μ Z γ EW from the original variance. The data statistical uncertainty has the largest impact on the measured signal strength, followed by the signal acceptance uncertainties. A small correlation is observed among the different sources of uncertainty. The signal uncertainties are divided into acceptance uncertainties for the signal events entering the fiducial volume and the shape uncertainties, which are the uncertainties in the shape of signal distributions within the fiducial volume. The signal acceptance uncertainties are assigned to the theoretical cross-section and not to the fiducial cross-section.
The measured fiducial cross-section is extracted by taking the product of the signal strength, μ Z γ EW , and the predicted cross-section times branching ratio to neutrinos in the fiducial volume defined in Sect. 6.2. The measurement Table 7 The contributions from different groups of systematic uncertainties to the ±1σ uncertainty bands of the μ Z γEW best-fit value and on B inv and B(H → γ γ d ) 95% CL limits. The evaluation is performed by fixing a given group of systematic uncertainties to their best-fit values and subtracting the new variance (σ 2 ) of the best-fit value or the limit from the nominal variance including all systematic uncertainties. Due to residual correlations between categories, the sum in quadrature of the systematic uncertainties can differ from the actual value. The uncertainty due to the finite number of data events ('Data stats.') is obtained by fixing all systematic uncertainties to their best-fit values. The sum of all systematic uncertainties is estimated by subtracting the statistical variance component from the total variance. The experimental uncertainties and the uncertainty related to the size of MC simulated samples ('MC stats.') are treated as separate categories. The V γ + jets theory entry includes the theoretical uncertainties in strong Z γ + jets, EW W γ + jets and strong W γ + jets production for μ Z γEW ; however, for B inv and B(H → γ γ d ), it also includes those from EW Z γ + jets. For the last two columns the impact of systematic uncertainties is computed from a fit to data with B inv = 0 or B(H → γ γ d ) = 0 for each respective column and SM prediction agree within the measurement uncertainties. The measured fiducial cross-section is σ fid. Z (→νν)γ EW = 1.31 ± 0.20(stat) ± 0.20(syst) fb, which includes the contribution from the interference term with the strong production of Z γ + jets. The interference computed through MadGraph is 2% in the fiducial volume and is treated as an uncertainty in the EW Z γ + jets cross-section. The theoretical MadGraph cross-section including the 0.3% NLO QCD K -factor correction from VBFNLO is 1.27 ± 0.01(stat)±0.17(LO QCD MadGraph scale)±0.03(pdf) fb = 1.27±0.17 fb. The jet-veto is not part of the fiducial phasespace definition; the loss in efficiency in simulation for this veto is 5%. Rev.Cen. CR, and Fake-e CR are included in the likelihood function definition to provide constraints on the background contribution to the SR. The signal contribution in the Z γ Rev.Cen. CR is taken into account in the statistical analysis. In the likelihood function definition, all the normalization factors β from Eq. (3) are fixed to one except the ones corresponding to the W γ + jets and fake-e background contributions. The μ Z γ EW normalization factor is fixed to one because the shape of this contribution and the one of the H → inv. signal are so similar that leaving the former unconstrained would largely affect the search sensitivity. The result of the EW Z γ + jets cross-section measurement described in Sect. 8.1 further supports this assumption. The β Z γ strong normalization factor is fixed to one because the CR yields are not large enough to reduce the theoretical uncertainties. In addition, the observed normalization in the EW Z γ + jets cross-section measurement is consistent with unity in Table 6. The likelihood is  The H → inv. signal is scaled to a B inv of 37%. The post-fit uncertainties include statistical, experimental, and theoretical contributions. The lower panel shows the ratio of data to the post-fit yields and also compares the pre-fit and post-fit background predictions Table 8 Data yields and background predictions, after the fit to 139 fb −1 of data with the B inv signal normalization set to zero, for the four DNN output score bins of the SR and the inclusive CRs. The uncertainties in the backgrounds are derived by the fit and include the effects of nuisance parameter constraints and the correlation of system-atic uncertainties; as a result, the uncertainties in the total background can be smaller than the sum of those in the single contributions. The predicted signal yields, assuming a 37% branching ratio for H → inv. are shown with their associated uncertainties. A dash '-' indicates less than 0.01 events

Process
Fake-e CR W Rev.Cen. CR for the B inv search, the third and fourth bins are merged into one bin covering 0.6-1.0 in DNN output score. The result of the maximumlikelihood fit to data in the 4 SR and 15 CR bins is shown in Fig. 4, with the best-fit model propagated in all the regions. The SR bin-by-bin yields and CR yields are shown in Table 8 for the background contribution, a benchmark H → inv. signal contribution, and recorded data. The fitted normalization of the sum of the EW and strong W γ + jets events relative to the SM prediction is 1.07 ± 0.18; the fake-electron normalization is not reported because no comparison with the SM predictions is possible. The level of agreement between the data and the prediction from background simulation is good overall, and is better after the fit, as shown in the lower panel of Fig. 4. No evidence of a new physics contribution is visible on top of the background prediction. The observed (expected) upper limit on B inv is 0.37 (0.34 +0. 15 −0.10 ) at 95% CL. The impact on the limit from different groups of uncertainties is shown in Table 7, and the results are evaluated in the same way as for observation of EW Z γ + jets, described in Sect. 8.1, except that B inv is fixed to zero.

Fit model and results for H → γ γ d search
In the search for a Higgs boson decaying into a γ γ d pair, the most powerful discriminating observable is the photon-E miss T transverse mass m T (γ , E miss T ), so this observable is used to search for this new physics signal. The events enter-ing the dedicated SR γ d (see Sect. 6) are separated into five m T bins, as described in Sect. 6.4. Because the relative contributions of H → γ γ d signal produced through ggF and VBF production vary with m jj , events are also separated into two m jj categories, those with m jj < 1 TeV and those with m jj ≥ 1 TeV. A total of ten bins in the SR as well as ten bins in each CR enter the likelihood function definition, which is equivalent to the one in Eq. (4) other than having a different number of bins in the SR and CRs and a different signal benchmark model for the interpretation in the context of the H → γ γ d search. The fixing of the β and μ Z γ EW normalization factors to unity in Eq. (4) is repeated to be consistent with Sect. 8.2, but allowing them to float in the fit does not change the results within rounding. The SR bin-by-bin yields and CR yields are shown in Table 9 for the background contribution and recorded data yields after a fit to the background-only contributions.
The result of the maximum-likelihood fit with the B(H → γ γ d ) signal normalization set to zero in the ten SR bins and four inclusive CRs is shown in Fig. 5. The CRs are shown inclusively to reduce the number of bins presented, but the same binning as the SR is used in the fit model. The data and predictions from background simulation agree within the reported uncertainties, apart from a small deficit in the data in the bins corresponding to the highest m T values. The pre-fit background predictions in the highest m T bins are pulled down by the fit to data, and uncertainties describing these differences, which increase in this high m T range, come from the interference between EW and strong Z (→ νν)γ + jets production as well as a comparison between MadGraph5_aMC@NLO and Sherpa simulation for the strong-production Z (→ νν)γ + jets background contributions. Overall, the level of agreement is better after the fit, as shown in the lower panel of Fig. 5; no evidence of a new physics contribution is visible on top of the background prediction. Figure 6 shows the distribution of m T (γ , E miss T ) in the inclusive SR γ d (no m jj split), and also shows the shape of a H → γ γ d signal for two different mass hypotheses compared with same post-fit background predictions and data as for Fig. 5. Thus the reasons for the change in the pre-fit predictions for the highest m T values are the same as described for Fig. 5. Good agreement between data and the background expectations in the m T (γ , E miss T ) distribution is observed also in the CRs considered in the statistical analysis.
The statistical analysis sets an observed (expected) upper limit on B(H → γ γ d ) of 0.018 (0.017 +0.007 −0.005 ) at 95% CL when considering both the VBF and ggF Higgs boson production mechanisms at a Higgs boson mass of 125 GeV. If considering a BSM scalar boson with a mass of 125 GeV produced through VBF, the observed (expected) upper limit on the cross-section times branching ratio is 0.064 pb (0.064 +0.030 −0.019 pb) at 95%. For such a BSM scalar boson pro-duced through ggF, the observed (expected) upper limit on the cross-section times branching ratio is 10.2 pb (7.3 +3.4 −1.9 pb) at 95%, which shows that the sensitivity is dominated by the VBF production mode.
The 95% CL limit on σ VBF × B(H → γ γ d ) has also been calculated for a VBF-produced Higgs boson with different mass hypotheses in the narrow width approximation (NWA), ranging from 60 GeV to 2 TeV, as shown in Fig. 7. The crosssection for a VBF-produced Higgs boson decreases rapidly with increasing boson mass, leading to smaller signal yields in the SR. The signal corresponding to a high-mass scalar mediator peaks towards high values of m T (γ , E miss T ), where the smaller background leads to good sensitivity despite the small expected signal.
The impact of various sources of uncertainty on the B(H → γ γ d ) upper limit is shown in Table 7, evaluated in the same way as for the H → inv. search, described in Sect. 8.2. The statistical uncertainty of the yields of data events in SR γ d has the largest impact on the limit determination. A negligible correlation is observed among the nuisance parameters corresponding to the different sources of uncertainty.

Conclusion
Data collected from 139 fb −1 of 13 TeV proton-proton collisions by the ATLAS experiment during the Run 2 of the LHC are scrutinized in a VBF-favoured signature of two forward jets, E miss T , and a photon, to provide constraints on several SM and BSM processes. The observation of SM EW Z γ + jets production is reported with an observed (expected) significance of 5.2σ (5.1σ ). The fitted normalization for the EW Z γ + jets process relative to the SM prediction is μ Z γ EW = 1.03 ± 0.25, corresponding to a measured crosssection of 1.31 ± 0.29 fb in the considered fiducial volume. A search for Higgs bosons decaying solely into invisible particles is performed in the same final-state signature. Because no significant excess is observed, 95% CL upper limits of 0.37 (0.34 +0. 15 −0.10 ) are set on the observed (expected) branching ratio to invisible particles. A search for Higgs bosons decaying into a photon and a dark photon is also performed, and the results exclude at 95% CL cross-section times branching ratio values ranging from 0.15 pb for a scalar mediator with a mass of 60 GeV to 3 fb for a scalar mediator with a mass of 2 TeV. For a Higgs boson mass of 125 GeV, the observed (expected) 95% CL upper limit on the Higgs boson branching ratio to γ γ d is 0.018 (0.017 +0.007 −0.005 ), the most stringent to date. binning, and for the inclusive CRs. The uncertainties in the backgrounds are derived by the fit and include the effects of nuisance parameter constraints and the correlation of systematic uncertainties; as a result, the uncertainties in the total background can be smaller than the sum of those in the single contributions. The predicted signal yields, assuming the SM production cross-section for a 125 GeV Higgs boson and A H → γ γ d signal is shown for two different mass hypotheses, 125 GeV and 500 GeV, and scaled to a B(H → γ γ d ) of 2% and 1%, respectively. The lower panel shows the ratio of data to the sum of all the background contributions, a comparison with the pre-fit background prediction, and the signal-to-background ratio shifted by 1.0 (to share the same vertical axis). Events with m T (γ , E miss T ) larger than the rightmost bin boundary are added to that bin Fig. 7 The 95% CL upper limit on the Higgs boson production cross-section times branching ratio to γ γ d for different VBF-produced scalar-mediator-mass hypotheses in the NWA. The theoretically predicted cross-section of a Higgs boson produced via VBF and with the B(H → γ γ d ) = 5% is superimposed on the ±1σ and ±2σ NNLO QCD+NLO EW uncertainty bands of the expected production cross-section limit Acknowledgements 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, Argentina; YerPhI, Armenia; ARC, Australia; Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors' comment: All ATLAS scientific output is published in journals, and preliminary results are made available in Conference Notes. All are openly available, without restriction on use by external parties beyond copyright law and the standard conditions agreed by CERN. Data associated with journal publications are also made available: tables and data from plots (e.g. cross section values, likelihood profiles, selection efficiencies, cross section limits, ...) are stored in appropriate repositories such as HEP-DATA (http://hepdata.cedar.ac.uk/). ATLAS also strives to make additional material related to the paper available that allows a reinterpretation of the data in the context of new theoretical models. For example, an extended encapsulation of the analysis is often provided for measurements in the framework of RIVET (http://rivet.hepforge.

org/).]
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permit-ted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3