Search for invisible Higgs boson decays in vector boson fusion at $\sqrt{s} = 13$ TeV with the ATLAS detector

We report a search for Higgs bosons that are produced via vector boson fusion and subsequently decay into invisible particles. The experimental signature is two energetic jets with $\mathcal{O}(1) $ TeV invariant mass and $\mathcal{O}(100) $ GeV missing transverse momentum. The analysis uses 36.1 fb$^{-1}$ of $pp$ collision data at 13 TeV recorded by the ATLAS detector at the LHC. In the signal region the 2252 observed events are consistent with the background estimation. Assuming a 125 GeV particle with Standard Model cross sections, the upper limit on the branching fraction of the Higgs boson decay into invisible particles is 0.37 at 95% CL where 0.28 was expected. This limit is interpreted in Higgs portal models to set limits on the WIMP-nucleon scattering cross section. We also consider invisible decays of additional scalar bosons with masses up to 3 TeV for which the upper limits on the cross section times branching fraction are in the range of $0.3-2.0$ pb.

The data sample corresponds to an integrated luminosity of 36.1 fb −1 of pp collisions at √ s = 13 TeV recorded by the ATLAS detector at the LHC in 2015 and 2016. The experimental signature of the signal VBF production process is a pair of energetic quark jets with a wide gap in pseudorapidity (η) corresponding to the O(1) TeV value of the invariant mass (m j j ) of the two highest-p T jets in the event. 1 The signature for the decay process is the O(100) GeV value of the missing transverse momentum (E miss T ) that corresponds to the Higgs boson p T . The VBF topology offers a powerful rejection of the strongly produced backgrounds of single vector boson plus two jets, and the multijet background. In this analysis, the Higgs production in the gluon fusion process is subdominant to VBF and is considered as part of the Higgs signal.
Direct searches for invisible Higgs decays look for an excess of events over Standard Model expectations. The absence of an excess is interpreted as an upper limit on the branching fraction of invisible decays (B inv ) assuming the Standard Model production cross section [20] of the 125 GeV Higgs boson. Other published results have targeted a variety of production mechanisms-gluon fusion, VBF, W or Z associated production [21-25]-to set upper limits on B inv . The best limits are from the statistical combination of search results for which ATLAS reported an observed (expected) limit of 0.25 (0.27) [26] and CMS 0. 24 (0.23) [27] at 95% confidence level (CL). For these combinations the single input with the highest expected sensitivity is VBF, the channel pursued here, for which ATLAS reported an observed limit of 0.28 [28] and CMS 0.43 [27], with an expected limit of 0.31 in both experiments.
Global fits to the measurements of visible decay channels of the Higgs boson place indirect constraints on the sum of the branching fraction to final states that are not detected using current reconstruction and analysis techniques plus the branching fraction to the invisible final states described above. For this sum, denoted by B , ATLAS reports an observed (expected) limit of 0.49 (0.48) [26] and CMS of 0. 57 (0.52) [29] with similar but not identical assumptions. A combination of ATLAS and CMS results gives 0.34 (0.39) [30]. As noted in Ref. [28], there is complementarity between the direct search for invisible Higgs decays and the indirect constraints from the global fits. A null result in the former and a non-zero result in the latter would point to undetected decays or incorrect model assumptions as the cause.
In this analysis, several changes and improvements are made with respect to the previous ATLAS paper on this topic [28]. The event selections are changed to retain a good sensitivity despite the higher pileup. The trigger and hadronic objects are defined considering the simultaneous pp collisions in the same and nearby bunch crossings (pileup) (Section 2). The leading backgrounds are simulated using state-of-the-art QCD predictions (Section 3). The analysis uses three bins in m j j to increase the signal sensitivity (Section 4). The Z νν estimation relies only on the Z ee and Z µµ control samples, and is not affected by theoretical uncertainties of the W-to-Z extrapolation (Section 5). The systematic uncertainties are evaluated separately for each m j j bin (Section 6). The search is repeated for other scalars with masses up to 3 TeV, which can easily be reinterpreted for models not considered in this Letter (Section 7). Several aspects of the analysis have not changed compared to the ATLAS Run-1 analysis-e.g., subdetector descriptions, transfer factor method, Higgs portal models-and details of these may be found in Ref. [28].
2 Detector, trigger, and analysis objects ATLAS is a multipurpose particle physics detector with a forward-backward symmetric cylindrical geometry consisting of a tracking system, electromagnetic and hadronic calorimeters, and a muon system [31].
The trigger to record the sample containing the potential signal process used a two-level E miss T algorithm with thresholds adjusted throughout the data-taking period to cope with varying levels of pileup [32, 33]. The -1 system used coarse-spatial-granularity analog sums of the calorimeter energy deposits to require E miss T > 50 GeV. The second-level system [34] used jets that are reconstructed from calibrated clusters of cell energies [35] and requires E miss T > 70-110 GeV depending on the luminosity and the pileup level. The trigger efficiency [36] for signal events is 98% for E miss T > 180 GeV when comparing the trigger selection with the offline E miss T definition that contains additional corrections.
The triggers to record the control samples for background studies used lepton and jet algorithms [37]. The samples with leptonic W and Z decays were collected with a single-electron or -muon trigger with p T > 24 GeV (26 GeV) and an isolation requirement in 2015 (towards the end of 2016). The sample of multijet events was collected using a set of low-threshold single-jet triggers with large prescale values to keep the event rate relatively low.
For each event, a vertex is reconstructed from two or more associated tracks (t) with p T > 400 MeV. If multiple vertices are present, the one with the largest t (p T,t ) 2 is taken as the hard-scatter primary vertex.
Leptons ( = e, µ) are identified to help characterize events with leptonic final states from decays of vector bosons. Since the signal process contains no leptons, such events are used for the background estimation described in Section 5. Electrons (muons) must have p T > 7 GeV, |η| < 2.47 (2.5), and satisfy an isolation requirement. Electrons are reconstructed by matching clustered energy deposits in the electromagnetic calorimeter to tracks from the inner detector [38,39] and muons by matching inner detector and muon spectrometer tracks [40]. All leptons must originate from the primary vertex.
Jets are reconstructed from topological clusters in the calorimeters using the anti-k t algorithm [41] with a radius parameter R = 0.4. Jets must have p T > 20 GeV and |η| < 4.5. The subset of jets with p T < 60 GeV and |η| < 2.4 are jet vertex tagged ( ) [42] to suppress pileup effects, using tracking and vertexing, with 92% efficiency.
Cleaning requirements help suppress non-collision backgrounds [43]. Fake jets due to noisy cells are removed by requiring a good fit to the expected pulse shape for each constituent calorimeter cell. Fake jets induced by beam-halo interactions with the LHC collimators are removed by requirements on their energy distribution and the fraction of their constituent tracks that originate from the primary vertex.
In events with identified leptons, the leptons can also be reconstructed as jets, which can potentially lead to double counting of objects. The lepton-jet overlap in ∆R distance2 is resolved sequentially as follows. If an electron is near a jet with ∆R < 0.2, the jet is removed to avoid the double counting of electron energy deposits. If a remaining jet is near an electron with 0.2 ≤ ∆R < 0.4, the electron is removed. If a muon is near a jet with ∆R < 0.4 and the jet is associated with at least (less than) three charged tracks with p T > 500 MeV, the muon (jet) is removed.
The E miss T variable is the magnitude of the negative vector sum of the transverse momenta, − i ì p T,i , where i represents both the "hard objects" and the "soft term." The hard objects consist of leptons and jets, which are individually reconstructed and calibrated; the list excludes pileup jets, which are removed by a requirement. The soft term is formed from inner detector tracks not associated with the hard objects, but matched to the primary vertex. In the search region, E miss T is mostly from the recoil against the dijet system.

The
procedure is intended to remove pileup jets, but can cause large fake E miss T if it removes a high-p T jet from the hard scatter, e.g., a jet from a p T -balanced three-jet event. In order to reduce this, a correlated quantity H miss T -defined as | j ì p T, j |, where j represents all jets without the requirement-is required to be H miss T > 150 GeV. In the three-jet example, H miss T would be near zero.
The E miss T significance (S ) is used only in events with an electron and is defined as where the p T quantities are for leading jet ( j 1 ), subleading jet ( j 2 ), and electron, respectively. The use of this quantity to reduce the contamination from jets misidentified as electrons is discussed in Section 5.

Event simulation
Monte Carlo simulation (MC) consists of an event generation followed by detector simulation [44]  , were subsequently reweighted to reproduce the pileup distribution in data. In general, simulated events were corrected for the trigger efficiency, the jet energy, and the lepton selections using dedicated data samples. Simulated events were corrected for the small differences between data and MC in the trigger, the lepton identification efficiency, and the jet energy scale and resolution.
For the signal process, the VBF cross sections were calculated at next-to-leading order (NLO) in QCD and the events were generated using -2 [49]; NLO electroweak corrections were applied using [50]. The generated events were interfaced with 8 [46] for hadronization and showering, using the tune [51] and the 3.0 NNLO PDF set [52]. The gluon fusion cross section was calculated at NNLO in QCD and events were generated using - [53] with the 4 15 PDF set [54] interfaced to a fast detector simulation [55]. The showering simulation followed the same procedure as for the VBF sample. For both the VBF and gluon fusion events, the H → Z Z * → 4ν process is included in the sample as invisible Higgs decays. Additional scalars with masses up to 3 TeV were simulated as described above for VBF signal process, assuming a full width of 4 MeV.
The W and Z events were generated using 2.2.1 [56] with [57] and [58] matrixelement generators, and merged with parton shower [59] using the + @ prescription [60]. The 3.0 NNLO PDF set was used. In terms of the order of the various processes, the strong production was calculated at NLO for up to two jets and leading order (LO) for the third and fourth jets. The electroweak production was calculated at LO for the second and third jets. The levels of the interference between electroweak and strong processes were computed with 5_a @ [61]. The interference on the total expected background is only 0.1% and thus neglected.
Other potential background processes involve top quarks, dibosons, and multijets. Top quarks and dibosons were generated with interfaced with and [62], which simulate the heavy-flavor decays. The diboson backgrounds include electroweak-mediated processes. The multijet estimate does not directly use the MC.

Event selection
All events must have a primary vertex. The selection listed below divides the data sample into a signalenriched search region (SR) and background-enriched control regions. The control regions and the statistical fit are discussed in detail in Section 6. The rest of this section focuses on the SR and the prefit event yields.
For the SR, an event is required to have • no isolated electron or muon, • a leading jet with p T > 80 GeV, • a subleading jet with p T > 50 GeV, • no additional jets with p T > 25 GeV, The two jets are required to have the following properties: • be in opposite η hemispheres, η j 1 · η j 2 < 0, • m j j > 1 TeV. Table 1: Event yields in the signal region (SR) and control regions (CR) summed over lepton charge and flavor. The yields are the prefit values for m j j > 1 TeV. The observed data (N), the background estimate (B), and the signal (S for m H = 125 GeV with B inv = 1) are given. The B and S values for individual processes are rounded to a precision commensurate with the sampling uncertainty associated with the finite MC sample size. For all processes the fractions of electroweak production [ ] are given. "Other" is defined in the text. The SR includes background events containing a W or Z plus two jets, where the W decays into eν, µν, and τν, and the Z decays into two neutrinos. Table 1 gives the prefit SR yields in the first column. The VBF production process gives the biggest contribution (87%) to the signal sample (fixed as B inv = 1). The contribution from gluon fusion accompanied by parton radiation is small (13%) and other production modes contribute negligibly. The fraction of VBF signal events that pass the signal region event selections, defined as acceptance times reconstruction efficiency, is 0.7%. For the backgrounds, both the strong production and the electroweak production contribute in the SR. The strong production processes contributes more than 70% of the backgrounds in all of the m j j bins.
As is discussed in Section 7, the signal significance is improved by considering three bins of m j j defined by boundaries at [1, 1.5, 2, -] TeV. The prefit S/B ratio (for B inv = 1) in these bins is approximately 0.3, 0.4, 0.8, respectively.

Control samples and statistical treatment
The main backgrounds in the SR are the W and Z processes and the minor backgrounds are the diboson, tt, and multijet processes. Accurate estimation of the W and Z processes is the biggest challenge of the analysis. Both estimations make use of control regions (CR) in the MC and the lepton-triggered data samples.
The W CR requires one identified lepton, but the selections are otherwise identical to those of the SR. The lepton-p T threshold is 30 GeV. The sample is divided into four subsamples depending on the lepton flavor and charge. The two W e ± ν subsamples are further subdivided by S < 4 √ GeV (> 4) to provide a subsample enriched (depleted) in fake electrons, where a jet is misidentified as an electron. In addition, a region enriched in non-prompt electrons was defined by requiring that the electron likelihood fail the tightest definition, while satisfying the looser definition. After subtracting the prompt W events, this sample was used to measure the ratio of events with S ≶ 4 √ GeV. The E miss T is calculated excluding all detector signals associated with leptons to mimic the quantity used in the SR. The kinematic bias in E miss T due to the S selection was found to be negligible.
The Z CR is based on the same selection criteria as the SR, but the lepton veto is replaced by the requirement of two same-flavor opposite-sign leptons with |m − m Z | < 25 GeV. The sample is studied separately for the two lepton flavors. The leading lepton-p T threshold is the same as above, and the subleading lepton-p T threshold is 7 GeV. The sample is divided by lepton flavor. The E miss T is calculated as for the W CR. Table 1 gives the prefit CR yields for m j j > 1 TeV for the W (Z) CR in the third (fourth) columns. These prefit yields are the inputs for the statistical fit described below. The samples are very pure, as the relative contribution of the W (Z) CR is 95% (99%) from W (Z) decays. The definitions of the parameters that are inputs to the fit are where the event yields are for the observed data (N) and the MC estimate of the background (B). The transfer factor α is the SR-to-CR ratio of the MC yields, and is a quantity useful for visualizing how the systematic uncertainties partially cancel out. The normalization β is the data-to-MC ratio in the CR, which is extracted from the fit.
Model testing uses a profile likelihood-ratio test statistic [63] in the CL -modified frequentist formalism [64]. A maximum-likelihood fit of the observed data and MC estimate for each bin sets an upper limit,3 using a one-sided CL, on B inv for the 125 GeV Higgs boson and on the product σ scalar · B inv for a scalar of different mass. The fit considers a total of 27 bins: three m j j bins for each of nine subsamples (one for the SR, six for the W CR, two for the Z CR). The prefit comparisons of data and MC are shown for all subsamples in Fig. 1.
Six normalization β parameters are extracted from the fit, one for each of the three m j j bins for the W and Z backgrounds. The β parameters extracted from the fit are consistent with unity within their 1 σ uncertainties. The β W (β Z ) parameters are extracted in a simultaneous fit of the six W CR (two Z CR) subsamples to the SR, one β for each m j j bin. In particular, the W eν subsamples are split into two bins of of S , one enriched in non-prompt electrons, and split by charge, since the non-prompt contribution is expected to be charge symmetric. The normalization of the fake component in the W eν subsamples with S ≶ 4 √ GeV are simultaneously determined, where the ratio between the fakes in the two regions is fixed. This ratio is determined by the fit using dedicated control region of electrons that satisfy a looser definition than is used in  1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 Fig. 2 for m j j and E miss T for the W and Z CR. The m j j (E miss T ) plot groups the backgrounds to show the dependence of the distribution shape on the production mechanism (final state).
The remaining processes-top quarks, dibosons, multijets-contribute negligibly to the SR (called "other" in Table 1). The first two are estimated with MC using nominal cross sections. The multijet contribution is very small, but it is a difficult process to estimate and a potentially dangerous background because those events that pass the E miss T selection are due to instrumental effects, such as the mismeasurement of the jet energy.
The billionfold-or-more reduction of multijets after the event selection makes it impractical to simulate, so a data-driven method based on a rebalance-and-smear strategy [66] is used. The assumption is that the E miss T is due to jet mismeasurement in the detector response to jets and neutrinos from heavy-flavor decays [67, 68]. Using the jet-triggered sample, the jet momenta are rebalanced by a kinematic fit, within their experimental uncertainties, to obtain the balanced value of the jets' p T . The rebalanced jets are smeared according to jet  Table 1, contribute a few events at low values of m j j and E miss T , and are omitted. The last bin in each plot contains the overflow.
response templates, which are obtained from MC and validated with dijet data. The procedure is validated in a ∆φ j j -sideband validation region (VR) with 95% purity. This VR is defined by 1.8 < | ∆φ j j | < 2.7 and the loosening of the other requirements (| ∆η j j | > 3, m j j > 0.6 TeV, and allow a third leading jet with 25 < p T < 50 GeV, but no other jets with p T > 25 GeV). The comparison of the predictions and the data in the VR shows good agreement (Fig. 3). The multijet component is obtained using the data-driven method with the associated systematic uncertainty bands, while the non-multijet components are obtained using MC.

Uncertainties
Experimental sources of uncertainty are due mainly to the jet energy scale and resolution [69], E miss T soft term [70], and lepton measurements [39,40]. In order to reduce fluctuations due to limited MC sample size, the uncertainties in number of expected events for the variations of jet energy scale and resolution for the strong and electroweak background samples are averaged. This is motivated by the similarities of the kinematics and the detector effects for the two production processes for each m j j bin. For the lepton measurements, the impact of identification in the W CR is negligible, but the veto in the SR affects the W background there. Other sources, such as the pileup distribution and luminosity [71,72], have a relatively small impact.
Theoretical sources of uncertainty are due mainly to scale choices in fixed-order matrix-element calculations. For the background processes, QCD scales are varied for the resummation scale (resum.), renormalization scale (renorm.), factorization scale (fact.), and matching scale. The first three scales in the listtechnically called q 2 , µ R , µ F , respectively-are varied by a factor of two. For the matching scale between the matrix element and the parton shower [56], the central value and the considered variations are 20 +10 −5 GeV. The higher-order electroweak corrections to the strongly produced W or Z are found to be negligible.
The effects of the theoretical variations are evaluated with a sample of generated MC events prior to reconstruction, which is larger than the reconstructed sample. Moreover, in order to reduce fluctuations due to limited MC statistics, the effect of the resummation and variations as a function of m j j are  determined by a linear fit, using m j j values below the selection for the SR and a sample with loosened selection on ∆η j j and ∆φ j j .
For both signal and background, the effects of the choice of a parton distribution function (PDF) set have a relatively small impact. The variations are considered using an ensemble of PDFs within the set [52] and the standard deviation of the distribution is taken as the uncertainty.
For the signal process, the effect of the scale uncertainty on the third-jet veto for the gluon fusion plus two-jet contribution is evaluated using the jet-bin method [73].
Statistical uncertainties are due to the data and MC sample sizes.
Systematic uncertainties are assumed to be either fully correlated or uncorrelated. The uncertainties from the following sources in each independent m j j bin are correlated between the SR and CR: QCD scales, PDF, and lepton measurements. The theoretical uncertainties due to QCD scales are uncorrelated between the following pairs: signal vs. background, electroweak vs. strong production, and W vs. Z production.
The sources of uncertainty are grouped into the three main categories given above ( Table 2). The impact of each source is measured in two ways: (1) on the 95% CL upper limit on B inv and (2) on the event yields and α transfer factors. Impact (1) assesses the percentage improvement of the B inv limit should that source of uncertainty be "removed" by fixing the associated parameter to its best-fit value. Impact (2) demonstrates that the systematic uncertainties in the individual yields partially cancel out for many of the theoretical sources. However, for many of the experimental sources the cancellation is not achieved due to limited MC statistics of the varied samples. For example, the effects of changing the renormalization and factorization scales change the MC yield in the Z SR (B Z ) and the Z CR (B Z ) by about 20%, but the α Z transfer factor changes by only 1%. In Table 2, only the 1 < m j j < 1.5 TeV yields are shown for the purpose of illustrating the partial cancellation. For the sources contributing the largest uncertainties, the α Z and the α W variations in the three m j j bins are shown graphically in Fig. 4. Table 2: Sources of uncertainty. The first set shows ∆, the relative improvement of the 95% CL upper limit on B inv when the source of uncertainty is "removed" by fixing it to its best-fit value. Combined rows are not simple sums of the rows above because of the ∆ metric; the symbols ( †, ‡, ) are parenthetically defined in the table. The column labeled "visual" shows bars whose lengths from the center tick are proportional to ∆. The second set of columns shows the effect on the yields and the α transfer factors; values in the 1 < m j j < 1.5 TeV bin are shown. The yields are for the signal process in the SR (S), Z MC in the SR (B Z ), and Z MC in the CR (B Z ). The α Z is given to demonstrate the reduction in the uncertainty in the ratio B Z /B Z . The individual yields for the W are not shown because the cancellation effects are similar to the Z counterparts. The abbreviations for the theoretical sources are described in the text. The '-' indicates that the quantity is not applicable. The penultimate (last) row shows the summary impact of removing the systematic uncertainties due to the experimental and theoretical sources (as well as statistical uncertainties of the MC samples).

Source
B inv improve.
[%] Yields, α changes (%) using all m j j bins in 1<m j j <1.5 TeV The combination of uncertainties from various sources shows that the dominant category has a systematic origin (penultimate row of Table 2). The lack of MC statistical precision for background processes with m j j > 2 TeV has the largest impact on B inv . We note that the ∆ values are percent improvements of B inv , so they do not add in quadrature or in any such standard statistical combinations.

Result and interpretations
The 2252 observed events in the SR are divided among the three m j j bins defined previously: 952, 667, and 633 events. These values are consistent with the background-only postfit yields of the sum of the background processes of 2100 events, which are divided among the three m j j bins: 850 ± 113, 660 ± 90, and 590 ± 812, respectively. The uncertainty represents the combined effect due to experimental and theoretical systematic uncertainties (MC sample size). These postfit values are also consistent with the prefit values. The expected signal yields (for B inv = 1 for VBF and gluon fusion) are 300, 310, and 460, respectively, and the last m j j bin, with S/B ≈ 0.8, has the highest sensitivity.
The postfit SR event distributions of m j j and E miss T are shown in Fig. 5, and a good agreement between the data and the expected backgrounds is observed. Figure 5(a) also shows that the S/B ratio rises with increasing m j j values, which motivates our division of the SR into multiple bins. The total electroweak contribution in the SR is relatively small at O(10%) ( Table 1), but the much flatter distribution of m j j makes it an important contribution to the final result. As noted in Section 5, the background estimation is done independently for each m j j bin to reduce the dependence on m j j modeling.
The fit, assuming the 125 GeV Higgs boson, gives the observed (expected) upper limit on B inv of 0.37 (0.28 +0. 11 −0.08 ) at 95% CL, and 0.32 (0.23 +0.11 −0.10 ) at 90% CL, where the uncertainties placed on the expected limit represent the 1σ variations. With this result, connections to dark matter can be made in the context of Higgs portal models [74]. In particular, relations between Higgs boson and scalar and Majorana fermion [11,75,76] allows the translation of the results into the -nucleon scattering cross section (σ -nucleon ).
The overlay of the interpretation of this result with the limits from some of the direct detection experiments [77][78][79] shows the complementarity in coverage (Fig. 6(a)). For the scalar interpretation cross sections are excluded at values ranging from O(10 −42 ) to O(10 −45 ) cm 2 and for the Majorana fermion interpretation the exclusion range is from O(10 −45 ) to O(10 −46 ) cm 2 , depending on the mass. The uncertainty band in the plot uses an updated computation of the nucleon form factors [80].
The correlation between B inv and σ -nucleon is presented in the effective field theory framework assuming that the new-physics scale is O(1) TeV [28], well above the scale probed at the Higgs boson mass. Adding a renormalizable mechanism for generating the fermion masses could modify the above-mentioned correlation [81].
In place of the 125 GeV Higgs boson, the same selection is applied to additional scalars with masses (m scalar ) of up to 3 TeV assuming only VBF production. The fraction of VBF signal events that pass the signal region event selections corresponding to the acceptance times efficiency ranges from 3-5%. The limit on σ scalar · B inv as a function of m scalar is shown in Fig. 6(b). The 95% confidence level upper limits on the cross section times branching fraction are in the range of 0.3-2.0 pb.

Conclusions
A search for Higgs boson decays into invisible particles is presented using the 36.1 fb −1 of pp collision data taken at √ s = 13 TeV by the ATLAS detector at the LHC. The targeted signature is the VBF topology with two energetic jets with a wide gap in η and large E miss T . Assuming the Standard Model cross section for the 125 GeV Higgs boson, an upper limit is set on B inv at 0.37 at 95% CL. This result is interpreted using Higgs portal models to exclude regions in the σ -nucleon vs.