Search for nonresonant pair production of Higgs bosons in the
$b\bar{b}b\bar{b}$ final state in $pp$ collisions at $\sqrt{s}=13$ TeV with
the ATLAS detector A search for nonresonant Higgs boson pair production in the
$b\bar{b}b\bar{b}$ final state is presented. The analysis uses 126 fb-1 of
$pp$ collision data at $\sqrt{s}={13}$ TeV collected with the ATLAS detector
at the Large Hadron Collider, and targets both the gluon–gluon fusion and
vector-boson fusion production modes. No evidence of the signal is found and
the observed (expected) upper limit on the cross-section for nonresonant Higgs
boson pair production is determined to be 5.4 (8.1) times the Standard Model
predicted cross-section at 95% confidence level. Constraints are placed on
modifiers to the $HHH$ and $HHVV$ couplings. The observed (expected) $2\sigma$
constraints on the $HHH$ coupling modifier, $\kappa_{\lambda}$, are determined
to be $[-3.5,11.3]$ ($[-5.4,11.4]$), while the corresponding constraints for
the $HHVV$ coupling modifier, $\kappa_{2V}$, are $[-0.0,2.1]$ ($[-0.1,2.1]$).
In addition, constraints on relevant coefficients are derived in the context
of the Standard Model effective field theory and Higgs effective field theory,
and upper limits on the $HH$ production cross-section are placed in seven
Higgs effective field theory benchmark scenarios. CERN-EP-2022-235 Phys. Rev.
D.

##  1 Introduction

The discovery of the $125\text{\,}\mathrm{GeV}$ Higgs boson ($H$) [1, 2, 3, 4]
at the Large Hadron Collider (LHC) has prompted a broad research program to
investigate its properties and compare the measurements with the Standard
Model (SM) predictions. Of particular interest is the search for nonresonant
Higgs boson pair, known as di-Higgs
($\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$),
production. This process has a strong dependence on the Higgs self-coupling,
which is a key ingredient of the electroweak symmetry breaking mechanism and a
sensitive probe for physics beyond the SM (BSM physics) in various scenarios,
such as two-Higgs-doublet models [5], composite Higgs models [6], twin Higgs
models [7], and the minimal supersymmetric extension of the SM [8, 9]. The
Higgs self-coupling also plays a fundamental role in understanding the
stability of the universe [10].

The dominant SM
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production process is gluon–gluon fusion (ggF). Its cross-section, for a Higgs
boson mass $m_{H}=$125\text{\,}\mathrm{GeV}$$, calculated at next-to-next-to-
leading order (NNLO) including finite top-quark-mass effects [11], is
$31.05\text{\,}\mathrm{f}\mathrm{b}$ at a center-of-mass energy
$\sqrt{s}=$13\text{\,}\mathrm{TeV}$$. The two dominant leading-order Feynman
diagrams contributing to this process are shown in Figure 1, where Figure 1(a)
is commonly referred to as the box diagram and Figure 1(b) as the triangle
diagram. The triangle diagram introduces the dependence on the trilinear Higgs
self-coupling, $\lambda$, shown by the red vertex in Figure 1(b), which can be
expressed in terms of its modifier, $\kappa_{\lambda}$.111A coupling modifier,
$\kappa$, is defined as the ratio of the modified coupling to its SM value,
$\kappa=c/c^{\text{SM}}$. By definition, $\kappa=1$ denotes the value of the
coupling predicted by the SM. In the SM, these two diagrams interfere
destructively. As a result, the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production cross-section and kinematic properties depend critically on the
value of $\kappa_{\lambda}$.

((a))

((b))

Figure 1: The two leading-order gluon–gluon fusion di-Higgs production Feynman
diagrams: (a) the box diagram; (b) the triangle diagram.

The
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production process with the second-highest cross-section in the SM is vector-
boson fusion (VBF), with a calculated value of
$1.73\text{\,}\mathrm{f}\mathrm{b}$ at next-to-next-to-next-to-leading order
(N3LO) [12], for $m_{H}=$125\text{\,}\mathrm{GeV}$$ at
$\sqrt{s}=$13\text{\,}\mathrm{TeV}$$. Figure 2 illustrates the Feynman
diagrams involved in di-Higgs production via vector-boson fusion at leading
order (LO). The coupling modifiers $\kappa_{\lambda}$, $\kappa_{V}$ and
$\kappa_{2V}$ are respectively shown at the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$,
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$VV
and
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$VV
interaction vertices, where $V$ stands for the gauge vector bosons
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
or
$\mathit{{Z}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$.
In the SM, the divergences in the Figures 2(b) and 2(c) diagrams exactly
cancel out due to perturbative unitarity. As $\kappa_{V}$ and $\kappa_{2V}$
depart from their SM value of one, this canceling out no longer occurs,
introducing a linear dependence of the cross-section on the effective center-
of-mass energy of the incoming vector bosons [13]. Therefore, the Higgs bosons
produced in non-SM $\kappa_{V}$/$\kappa_{2V}$ scenarios are expected to be
more energetic and more central in the detector on average. This increase in
the energy of Higgs bosons with increasing deviation from the SM continues up
to the scale of some new physics, which is required to unitarize the total
amplitude.

((a))  ((b))

((c))

Figure 2: The three tree-level vector-boson fusion di-Higgs production Feynman
diagrams.

The analysis described in this paper targets the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
process in the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
final state, in both the ggF and VBF production modes, using the data
collected by ATLAS between 2016 and 2018, during Run 2 of the LHC. Assuming
the SM branching ratio of 58.2% for
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
[14, 15], about one third of di-Higgs events decay into
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$,
making it the most abundant di-Higgs final state. However, as this is a fully
hadronic final state, the analysis faces the challenge of large backgrounds,
which originate mostly from nonresonant QCD production of multiple heavy
($\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}/\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$)
quarks, as well as from light-quark-initiated jets misidentified as
originating from heavy quarks.

The results are interpreted in terms of constraints on the $\kappa_{\lambda}$
and $\kappa_{2V}$ coupling modifiers, assuming $\kappa_{V}=1$. The analysis
also provides one- and two-dimensional constraints on relevant couplings in
the SM effective field theory (SMEFT) [16, 17, 18] and Higgs effective field
theory (HEFT) [19, 20] frameworks. In the SMEFT framework, the effects of new
physics may be described with an effective Lagrangian:

$\mathcal{L}_{\text{SMEFT}}=\mathcal{L}_{\text{SM}}+\frac{1}{\Lambda^{2}}\sum_{k}c_{k}^{(6)}O_{k}^{(6)}\,,$

where $\mathcal{L}_{\text{SM}}$ represents the SM Lagrangian, $O_{k}$ are
higher-dimensional local operators, $c_{k}$ are the Wilson coefficients, and
$\Lambda$ is the mass scale of the new physics phenomena (set to
$1\text{\,}\mathrm{TeV}$ for this result). The analysis considers operators
$O_{k}$ in the Warsaw basis, which provides a complete set of operators
allowed by SM gauge symmetries at dimension six [21] (dimension-five operators
introduce lepton and baryon number violation, and are therefore ignored in
this result). The five operators relevant to the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
process and their coefficients,
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\square}$,
$c_{t\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
$c_{tG}$, and
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}G}$,
are listed in Table 1 [22]. The computation of amplitudes from the above
Lagrangian includes three terms: a pure SM term, a “quadratic” term of order
$(1/{\Lambda^{4}})$ including purely new physics, and a “linear” term of order
$(1/\Lambda^{2})$ accounting for the interference between the SM and new
physics. The SMEFT constraints calculated in this analysis include both the
linear and quadratic new physics terms.

Table 1: The five relevant SMEFT coefficients and their corresponding dimension-6 operators, as defined in the Warsaw basis [21, 22]. Wilson Coefficient | Operator
---|---
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$ | $(\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}^{{\dagger}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})^{3}$
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\square}$ | $(\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}^{{\dagger}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})\square(\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}^{{\dagger}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})$
$c_{\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$ | $(\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}^{{\dagger}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})(\bar{Q}\tilde{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}t)$
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}G}$ | $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}^{{\dagger}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}G^{A}_{\mu\nu}G_{A}^{\mu\nu}$
$c_{\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}G}$ | $(\bar{Q}\sigma^{\mu\nu}T^{A}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})\tilde{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}G^{A}_{\mu\nu}$

In the HEFT framework, new physics in the electroweak sector is described
through anomalous couplings of the Higgs boson. The organization of the HEFT
Lagrangian is guided by chiral perturbation theory [23], with the low-energy
dynamics of electroweak symmetry breaking described using a nonlinear
realization of the gauge symmetry group $SU(2)_{L}\times U(1)_{Y}$. One
advantage of the HEFT framework is that the anomalous single-Higgs-boson and
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
couplings are defined separately, allowing simplified
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
interpretations. In the HEFT Lagrangian, ggF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production is described at LO by five relevant operators and their associated
Wilson coefficients:
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
$c_{t\overline{t}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
$c_{gg\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
$c_{gg\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
and
$c_{t\overline{t}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$.
In this formalism,
$c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
is equivalent to $\kappa_{\lambda}$ and
$c_{t\overline{t}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
is equivalent to the modifier for the coupling between the Higgs boson and top
quark, $\kappa_{t}$, shown by the light blue vertex in Figure 1. Fixing
$c_{t\overline{t}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=c_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=1$
and
$c_{gg\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=c_{gg\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=c_{t\overline{t}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=0$
restores the SM. At next-to-leading order (NLO), seven HEFT benchmark models
(BM) [24] have been defined using cluster analysis [25] to probe a wide
variety of characteristic shapes of the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
spectrum resulting from different BSM scenarios. The values of the
coefficients used to define these scenarios are given in Table 2.

Table 2: The values of the HEFT Wilson coefficients in the SM and in seven BSM benchmark models, as defined in Ref. [24]. Benchmark Model | $c_{HHH}$ | $c_{ttH}$ | $c_{ggH}$ | $c_{ggHH}$ | $c_{ttHH}$
---|---|---|---|---|---
SM | 1 | 1 | 0 | 0 | 0
BM1 | $3.94$ | $0.94$ | $1/2$ | $1/3$ | $-1/3$
BM2 | $6.84$ | $0.61$ | $0.0$ | $-1/3$ | $1/3$
BM3 | $2.21$ | $1.05$ | $1/2$ | $1/2$ | $-1/3$
BM4 | $2.79$ | $0.61$ | $-1/2$ | $1/6$ | $1/3$
BM5 | $3.95$ | $1.17$ | $1/6$ | $-1/2$ | $-1/3$
BM6 | $5.68$ | $0.83$ | $-1/2$ | $1/3$ | $1/3$
BM7 | $-0.10$ | $0.94$ | $1/6$ | $-1/6$ | $1$

The ATLAS Collaboration has previously published search results for
nonresonant
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
production using $27\text{\,}{\mathrm{fb}}^{-1}$ of early Run 2 data [26], and
a dedicated search for VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production in $126\text{\,}{\mathrm{fb}}^{-1}$ of data collected between 2016
and 2018 [27]. The present analysis benefits from the use of the 2016–2018
data for both production channels and also takes advantage of improvements in
jet reconstruction and in the identification of jets arising from the decays
of
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-quarks
(“$b$-tagging”) achieved by the ATLAS Collaboration since the publication of
Ref. [26]. In addition, the analysis employs a fully data-driven technique for
the background estimation, which uses an artificial neural network to perform
a kinematic reweighting of data from an alternative phase space to model the
background in the region of interest. The CMS Collaboration has also published
results of a search for nonresonant
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
with its full Run 2 dataset [28], setting the observed (expected) upper limit
on the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
cross-section at 3.9 (7.8) times the SM predicted cross-section, and
restricting the allowed interval for $\kappa_{\lambda}$ to $[-2.3,9.4]$
($[-5.0,12.0]$), both at 95% confidence level (CL). A more recent CMS
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
publication [29], in which the analysis exploits topologies arising from
highly energetic Higgs boson decays into $b\bar{b}$, sets the observed
(expected) upper limit at 9.9 (5.1) times the SM cross-section expectation,
and restricts the allowed interval for $\kappa_{2V}$ to [0.62, 1.41] ([0.66,
1.37]), at 95% CL. Other searches for nonresonant $HH$ production were
performed by ATLAS and CMS in the $b\bar{b}\tau^{+}\tau^{-}$ [30, 31],
$b\bar{b}\gamma\gamma$ [32, 33], $b\bar{b}\ell^{+}\nu\ell^{-}\nu$ [34, 35]
decay channels, as well as by ATLAS in the $b\bar{b}qq\ell\nu$ [36],
$WW^{*}\gamma\gamma$ [37] and $WW^{*}WW^{*}$ [38] decay channels. Amongst
them, the most sensitive results to date from ATLAS come from the
$b\bar{b}\gamma\gamma$ analysis, which sets the observed (expected) 95% CL
upper limit on the SM nonresonant $HH$ cross-section at 4.2 (5.7) times the SM
expectation and restricts the corresponding $\kappa_{\lambda}$ interval to
$[-1.5,6.7]$ ($[-2.4,7.7]$). The most sensitive results to date from CMS come
from the combination of the $b\bar{b}ZZ$, multilepton, $b\bar{b}\gamma\gamma$,
$b\bar{b}\tau\tau$, and
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
analyses, which set the observed (expected) 95% CL upper limit on the SM
nonresonant $HH$ cross-section at 3.4 (2.5) times the SM expectation and
restricts the corresponding observed $\kappa_{\lambda}$ interval to
$[-1.24,6.49]$ [39].

This document is structured as follows. The ATLAS detector and the data and
simulated events used in the analysis are described in Sections 2 and 3,
respectively. Section 4 presents the reconstruction and identification of
physics objects in this analysis and Section 5 details the event selection and
categorization. The background modeling method is described in Section 6, the
systematic uncertainties are detailed in Section 7 and, finally, the results
are reported in Section 8 and the conclusion is given in Section 9.

##  2 ATLAS detector

The ATLAS detector [40] at the LHC covers nearly the entire solid angle around
the collision point.222ATLAS uses a right-handed coordinate system with its
origin at the nominal interaction point (IP) in the center of the detector and
the $z$-axis along the beam pipe. The $x$-axis points from the IP to the
center of the LHC ring, and the $y$-axis points upwards. Cylindrical
coordinates $(r,\phi)$ are used in the transverse plane, $\phi$ being the
azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of
the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. Angular distance is
measured in units of $\Delta R\equiv\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$.
It consists of an inner tracking detector surrounded by a thin superconducting
solenoid, electromagnetic and hadron calorimeters, and a muon spectrometer
incorporating three large superconducting air-core toroidal magnets.

The inner-detector (ID) system is immersed in a $2\text{\,}\mathrm{T}$ axial
magnetic field and provides charged-particle tracking in the range
$|\eta|<2.5$. The high-granularity silicon pixel detector covers the vertex
region and typically provides four space-point measurements per track, the
first hit normally being in the insertable B-layer installed before Run 2 [41,
42]. Following the pixel detector is the silicon microstrip tracker, which
usually provides eight measurements per track. These silicon detectors are
surrounded by the transition radiation tracker, which enables radially
extended track reconstruction up to $|\eta|=2.0$.

The calorimeter system covers the pseudorapidity range $|\eta|<4.9$. Within
$|\eta|<3.2$, electromagnetic calorimetry is provided by barrel and endcap
high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin
LAr presampler covering $|\eta|<1.8$ to correct for energy loss in material
upstream of the calorimeters. Hadron calorimetry is provided by the
steel/scintillator-tile calorimeter, segmented into three barrel structures
within $|\eta|<1.7$, and two copper/LAr hadron endcap calorimeters. The solid
angle coverage is completed with forward copper/LAr and tungsten/LAr
calorimeter modules optimized for electromagnetic and hadronic energy
measurements respectively.

The muon spectrometer (MS) comprises separate trigger and high-precision
tracking chambers measuring the deflection of muons in a magnetic field
generated by the superconducting air-core toroidal magnets. The field integral
of the toroids ranges between $2.0$ and
$6.0\text{\,}\mathrm{T}$$\cdot$$\text{\,}\mathrm{m}$ across most of the
detector. A set of precision chambers covers the region $|\eta|<2.7$ with
three layers of monitored drift tubes, complemented by cathode-strip chambers
in the forward region, where the background is highest. The muon trigger
system covers the range $|\eta|<2.4$ with resistive-plate chambers in the
barrel, and thin-gap chambers in the endcap regions.

Interesting events are selected by the first-level trigger system implemented
in custom hardware, followed by selections made by algorithms implemented in
software in the high-level trigger [43]. The first-level trigger accepts
events from the $40\text{\,}\mathrm{MHz}$ bunch crossings at a rate below
$100\text{\,}\mathrm{kHz}$, which the high-level trigger reduces in order to
record events to disk at about $1\text{\,}\mathrm{kHz}$.

An extensive software suite [44] is used in data simulation, in the
reconstruction and analysis of real and simulated data, in detector
operations, and in the trigger and data acquisition systems of the experiment.

##  3 Data and simulated samples

###  3.1 Data sample

This analysis is performed in LHC proton–proton
($\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$)
collision data at $\sqrt{s}=$13\text{\,}\mathrm{TeV}$$ collected between 2016
and 2018. Only data collected during stable beam conditions are used, with all
relevant detector systems functional [45], corresponding to an integrated
luminosity of $126\text{\,}{\mathrm{fb}}^{-1}$. During 2016 data taking, a
fraction of the data ($8.3\text{\,}{\mathrm{fb}}^{-1}$) was affected by an
inefficiency in the online primary vertex reconstruction, which reduced the
efficiency of the $b$-tagging algorithms in the trigger; those events were not
retained for further analysis, resulting in an integrated luminosity of
$24.6\text{\,}{\mathrm{fb}}^{-1}$ for the 2016 dataset. The integrated
luminosities of the 2017 and 2018 datasets are
$43.7\text{\,}{\mathrm{fb}}^{-1}$ and $57.7\text{\,}{\mathrm{fb}}^{-1}$,
respectively.

The analysis uses events that satisfy either of two types of trigger
signatures, each with different requirements on the number of jets and their
$b$-tagging status [46]. The jets used are reconstructed with the anti-$k_{t}$
algorithm [47, 48], with a radius parameter of $R=0.4$. The $b$-tagging is
performed at the trigger level with the MV2c20 algorithm in 2016 and the
MV2c10 algorithm in 2017 and 2018 [46], with a range of
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-jet
identification efficiency operating points from 40% to 70% (as calculated from
simulated
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
samples.) The first of the two trigger signatures used for selecting
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
events requires two $b$-jets plus one additional jet (“2$b$1j”), while the
second requires two $b$-jets plus two additional jets (“2$b$2j”). The minimum
transverse energy ($E_{\text{T}}$) requirement on the jets is
$35\text{\,}\mathrm{GeV}$ for all jets used in the 2$b$2j trigger. In the
2$b$1j trigger, the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets must have $E_{\text{T}}>$55\text{\,}\mathrm{GeV}$$, while the requirement
on the minimum $E_{\text{T}}$ of the additional jet is between 100 and
$150\text{\,}\mathrm{GeV}$, depending on the year of data taking.

###  3.2 Simulated samples

Monte Carlo (MC) simulation is used for the modeling of signal events, as well
as to produce event samples of background processes for cross-checks and
validation studies. The Higgs boson mass is set to $125\text{\,}\mathrm{GeV}$
in the simulation. All samples were processed by the ATLAS simulation
framework [49] and the detector response was simulated with Geant4 [50].

The ggF signal process was simulated using the Powheg Box v2 generator [51,
52, 53] at NLO, including finite top-quark-mass effects, using the PDF4LHC15
[54] parton distribution function (PDF) set. Parton showers and hadronization
were simulated with Pythia 8.244 [55] with the A14 set of tuned parameters
[56] and the NNPDF2.3lo PDF set [57]. The SM ggF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
cross-section was taken as $\sigma_{\text{ggF}}=$31.05\text{\,}\mathrm{fb}$$,
calculated at NNLO including finite top-quark-mass effects [11]. Signal
samples for the ggF process were generated explicitly for coupling modifier
values of $\kappa_{\lambda}=1$ and $10$. A reweighting method is used to
obtain a ggF signal sample at each $\kappa_{\lambda}$ value, as described in
Ref. [58]: scale factors are derived as a function of $\kappa_{\lambda}$ in
bins of the generator-level invariant mass of the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
system by performing a linear combination of generator-level samples at three
different $\kappa_{\lambda}$ values ($\kappa_{\lambda}=0,1,$ and $20$). The
$\kappa_{\lambda}=10$ ggF signal sample is used to validate the derived scale
factors; this generated sample and the signal sample obtained from the
reweighting method are found to agree within the statistical precision of the
simulated sample. Additional generator-level ggF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal samples without parton showering were produced with Powheg Box v2 for
the $\kappa_{\lambda}=0$ and $20$ coupling modifier configurations to provide
a basis for the $\kappa_{\lambda}$ reweighting, along with the SM ggF sample.
For the reweighted ggF signal, the NNLO cross-section as a function of
$\kappa_{\lambda}$ is taken from Ref. [11]. In order to assess parton
showering uncertainties, alternative ggF samples were generated using the
Powheg Box v2 generator at NLO with the PDF4LHC15 PDF set, interfaced to
Herwig 7.1.6 [59] for parton showering and hadronization using the Herwig
7.1-default set of tuned parameters [60] and MMHT2014lo PDF set [61].

To extract SMEFT coefficient constraints, parton-level ggF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
samples were generated with MadGraph5_aMC@NLO [62, 63, 64] with the SMEFT@NLO
model [65] for a variety of SMEFT coefficients. A finely spaced
multidimensional grid of signal samples was obtained using a LO-derived
reweighting procedure in the generator-level invariant mass of the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
system; this procedure is similar to that used to obtain $\kappa_{\lambda}$
variations for the ggF signal, as described above. To extract HEFT coefficient
constraints, a similar NLO-derived reweighting procedure was applied to the
simulated reconstruction-level ggF signal sample to produce a variety of HEFT
signal scenarios, including the seven benchmark scenarios defined in Section
1, following the prescription outlined in Refs. [66, 67]. Additional
$K$-factors were applied to the SMEFT samples; these $K$-factors were derived
using the ratio of the NLO cross-section to the LO cross-section at the
equivalent HEFT point, as obtained using the HEFT to SMEFT translation from
Ref. [24].333Variations in the $c_{tG}$ Wilson coefficient were neglected when
calculating $K$-factors because the corresponding chromomagnetic operator does
not appear at LO within HEFT.

The VBF signal process was simulated using MadGraph 2.7.3 [63] at LO with the
NNPDF3.0nlo PDF set [68], interfaced with Pythia 8.244 for parton showering
and hadronization using the A14 set of tuned parameters and NNPDF2.3lo PDF
set. Signal samples for the VBF process were generated explicitly for coupling
modifier values of $(\kappa_{\lambda},\kappa_{2V},\kappa_{V})=(1,1,1)$,
$(1,1.5,1)$, $(2,1,1)$, $(10,1,1)$, $(1,1,0.5)$, $(-5,1,0.5)$, $(0,1,1)$,
$(1,0,1)$, and $(1,3,1)$. A linear combination of the first six of the listed
samples is used to derive distributions for a finer granularity of
$\kappa_{2V}$ values, following a technique used previously to generate
$\kappa_{\lambda}$ distributions [69]. The specific basis of six samples
utilized is chosen to avoid large statistical uncertainties in the reweighted
signal samples resulting from sparsely populated areas of kinematic phase
space. The generated VBF signal samples not included in the linear combination
basis – $(\kappa_{\lambda},\kappa_{2V},\kappa_{V})=(0,1,1)$, $(1,0,1)$, and
$(1,3,1)$ – were used to validate the performance of the combination method.
These generated samples and the corresponding signal samples obtained from the
combination method were found to agree within the statistical precision of the
simulated samples. The cross-section for the VBF $HH$ process, evaluated at
N3LO in QCD, is $1.73\text{\,}\mathrm{fb}$ in the SM [70, 71, 12, 72]. For the
reweighted VBF signal points, the N3LO to LO cross-section ratio at the SM
value is calculated, and this factor is applied to the cross-sections at each
$\kappa_{\lambda}$, $\kappa_{2V}$, and $\kappa_{V}$ point. In order to assess
parton showering uncertainties, alternative LO samples were generated using
MadGraph 2.7.3 with the NNPDF3.0nlo PDF set, interfaced to Herwig 7.0.4 with
the Herwig 7.1-default set of tuned parameters and MMHT2014lo PDF set for
parton showering and hadronization.

Top-quark pair production
($\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$)
and multijet background processes were simulated in order to validate the
background modeling procedure. The
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
sample was simulated at NLO in $\alpha_{\text{s}}$ using Powheg Box v2 [73].
Parton showering, hadronization, and the underlying event were modeled using
Pythia 8.230. The matrix element calculation uses NNPDF3.0nlo as the PDF set,
while the parton shower and underlying-event modeling uses NNPDF2.3lo and the
A14 set of tuned parameters. The damping parameter $h_{\mathrm{damp}}$, which
effectively regulates radiation at high $p_{\text{T}}$, was set to 1.5 times
the top quark’s mass. The
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
simulation is normalized using the value of the inclusive cross-section
calculated with Top++ 2.0 [74, 75]. This accounts for NNLO corrections in
$\alpha_{\text{s}}$, including next-to-next-to-leading logarithmic (NNLL)
resummation of soft gluon terms. The multijet background samples were modeled
using Pythia 8.235. This simulates pure QCD 2-to-2 interactions at LO in
$\alpha_{\text{s}}$. Events were showered using the parton shower native to
Pythia, which includes radiation and splitting that can result in additional
jets. The A14 set of tuned parameters and the NNPDF2.3lo PDF set were used.

Other background processes, such as SM Higgs boson,
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
(in other final states) and electroweak diboson production, have been
estimated to give negligible contributions to the selected event yields and
are therefore not included.

The effect of multiple interactions in the same and neighboring bunch
crossings (pileup) was modeled by overlaying each simulated hard-scattering
event with inelastic
$\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events generated with Pythia 8.186 using the NNPDF2.3lo PDF set and the A3 set
of tuned parameters [76]. Additionally, for all
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal samples, heavy-flavor decays were modeled using EvtGen 1.7.0 [77].

##  4 Object reconstruction

Primary vertices from
$\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{p}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
interactions are reconstructed [78] using at least two charged-particle tracks
with transverse momentum ($p_{\text{T}}$) above $500\text{\,}\mathrm{MeV}$
measured with the ID. The vertex with the largest sum of squared track momenta
($\sum p_{\text{T}}^{2})$ is taken as the hard-scatter primary vertex.

Hadronic jets are reconstructed using the anti-$k_{t}$ algorithm with radius
parameter $R=0.4$. The jet clustering uses particle-flow objects as inputs
[79]. Particle-flow objects are charged-particle tracks matched to the hard-
scatter vertex and calorimeter energy clusters after applying an energy
subtraction algorithm that removes the calorimeter deposits associated with
good-quality tracks from any vertex. The tracking information helps to improve
the energy resolution of the calorimeter clusters and reduce the impact from
pileup. The momenta of reconstructed jets are calibrated in a multistep
procedure [80]. Jets with $p_{\text{T}}<$60\text{\,}\mathrm{GeV}$$ and
$|\eta|<2.4$ must also satisfy a requirement based on the output of the
multivariate “jet vertex tagger” (JVT) algorithm [81], which is used to
identify and reject jets in which much of the energy originates from pileup
interactions. Correction factors are applied to the simulated events to
compensate for differences between the JVT efficiencies in data and
simulation. In the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
analysis, jets are discarded if they fail the “Tight” JVT working point,
corresponding to an average efficiency of 96% for jets from the hard-scatter
vertex.

Jets with radius parameter $R=0.4$ are also reconstructed from topological
clusters of energy deposits in the calorimeter [82] and calibrated in the same
way as the jets reconstructed from particle-flow objects. These jets are used
exclusively for the purpose of applying quality criteria to identify events
which are consistent with noise in the calorimeter or noncollision background
[83]. Events containing at least one such jet with
$p_{\text{T}}>$20\text{\,}\mathrm{GeV}$$, satisfying the JVT requirement, but
not these quality criteria, are rejected.

The identification of jets originating from $b$-quarks is performed by the
DL1r algorithm [84], which is applied to all jets with $|\eta|<2.5$. DL1r is
based on a multivariate classification technique combining information from
the impact parameters of ID tracks, the presence of displaced secondary
vertices, and the reconstructed flight paths of $b$\- and $c$-hadrons inside
the jet. The DL1r working point used in the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
analysis is the one that gives 77% efficiency for jets associated with true
$b$-hadrons in simulated
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
events. At this working point, the light-jet (charm-jet) rejection measured in
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
simulation is about a factor of 130 (4.9). The calibration of the DL1r
algorithm is performed separately for each jet type [85, 86] and correction
factors are derived and applied to the simulated samples to compensate for
differences between the $b$-tagging efficiencies in data and simulation.

Muons are reconstructed by matching ID tracks with either MS tracks or aligned
individual hits in the MS and performing a combined track fit [87]. They are
required to have $p_{\text{T}}>$4\text{\,}\mathrm{GeV}$$ and $|\eta|<2.5$, and
to satisfy “Medium” identification criteria based on track-quality variables.
Muons are used only to apply energy corrections to jets.

A momentum correction is applied to $b$-tagged jets to account for energy lost
to soft out-of-cone radiation and to muons and neutrinos in semileptonic
$b$-hadron decays. This correction follows the procedure used in Ref. [88] and
consists of two steps. First, a search is performed for muons located near the
jet which fall within a cone of variable size $\Delta
R(\mu,\text{jet})<\min\left(0.4,0.04+10/p_{\text{T}}^{\mu}\ \text{GeV}\right)$
around the jet axis. If a muon is found, its four-momentum is added to that of
the jet, and the energy deposited in the calorimeter by the muon is subtracted
from the jet to avoid double counting; this is computed according to the
description in Ref. [89]. In the second step, a global scale factor is applied
to each $b$-tagged jet according to its $p_{\text{T}}$ and whether or not it
has a muon associated with it. These scale factors are derived from
simulation.

##  5 Analysis selection and categorization

The analysis utilizes a set of criteria to select
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
candidate events, including dedicated requirements to separate events into
orthogonal ggF and VBF signal regions. “Forward” and “central” jets are used
with the following selection criteria:

  * •

central jets: $|\eta|<2.5$ and $p_{\text{T}}>$40\text{\,}\mathrm{GeV}$$;

  * •

forward jets: $2.5<|\eta|<4.5$ and $p_{\text{T}}>$30\text{\,}\mathrm{GeV}$$.

An initial “preselection” is applied to all events, which requires at least
four central jets with $p_{\text{T}}>$40\text{\,}\mathrm{GeV}$$, at least two
of which are $b$-tagged. As described in Section 3, the events considered in
this analysis are selected online through the 2$b$2j or 2$b$1j trigger
signatures. In order to simplify the modeling of trigger efficiencies, a
further selection is applied using offline kinematic quantities. Events are
selected if they have a leading444In this document, terms like “leading”,
“subleading” etc for physics objects refer to the ordering of these objects in
decreasing $p_{\text{T}}$. jet with $p_{\text{T}}>$170\text{\,}\mathrm{GeV}$$,
a third leading jet with $p_{\text{T}}>$70\text{\,}\mathrm{GeV}$$, and pass
the 2$b$1j trigger, or if they fail either of the two jet-$p_{\text{T}}$
requirements and pass the 2$b$2j trigger. This selection step leads to about a
10% loss of signal efficiency, but enables the reliable calculation of
simulation-to-data correction factors for estimating the trigger efficiency in
the remaining
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
signal events, depending on which of the above two trigger classes they belong
to.

Events passing the above preselection are required to contain at least four
central jets passing the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagging
requirement outlined in Section 4. The four highest-$p_{\text{T}}$
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets are chosen to reconstruct the decays of the two Higgs bosons. In about
$75\text{\,}$ of simulated signal events reaching this selection stage, these
four jets can be matched one-to-one (within $\Delta R<0.3$) to the four
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-quarks
from the decays of the Higgs bosons. In signal events where this matching
fails, one of the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-quarks
from the Higgs boson decays typically produces a jet that is outside the
analysis acceptance.

From the four selected
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets, there are three possible combinatorial pairings to form the two Higgs
boson candidates. Of those three configurations, the analysis selects the one
in which the higher-$p_{\text{T}}$ jet pair has the smallest $\Delta R$
separation. In the simulated samples with SM coupling values, for which the
analysis was mainly optimized, this method gives the correct pairing in around
90% of those signal events in which the four
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets are correctly matched to the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-quarks
from the decays of the Higgs bosons. While the pairing accuracy drops for
values of the coupling modifiers $\kappa_{\lambda}$ and $\kappa_{2V}$ that
result in softer $p_{\text{T}}$ spectra for the produced Higgs bosons, this
pairing method leads to a smoothly varying distribution of the expected
background in the plane of the invariant masses of the two Higgs boson
candidates, which facilitates the data-driven background estimation described
in Section 6.

Events are then subjected to additional selections designed to separate out
those consistent with the VBF production mode. For this, events must contain
at least two additional jets, central or forward;
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets are excluded. The two jets forming the pair with the largest invariant
mass ($m_{\text{jj}}$) are chosen as the “VBF jets”. The VBF jet pair is
required to satisfy $m_{\text{jj}}>$1\text{\,}\mathrm{TeV}$$, and the
pseudorapidity separation between the two jets, $|\Delta\eta_{\text{jj}}|$,
must satisfy $|\Delta\eta_{\text{jj}}|>3$. Lastly, the transverse component of
the momentum vector sum of the two VBF jets and the four jets forming the
Higgs boson candidates is required to be less than $65\text{\,}\mathrm{GeV}$.
Events satisfying the above criteria enter the VBF signal region, while those
failing to satisfy any of these criteria are considered further in the ggF
signal region.

Events satisfying either the ggF or VBF selections are required to satisfy
additional selection criteria designed to reduce the background and improve
the analysis sensitivity. In order to suppress the
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
background, a top-veto discriminant
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
is defined as:

$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=\text{min}\left[\sqrt{\left(\frac{m_{\text{jj}}-m_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}}{0.1m_{\text{jj}}}\right)^{2}+\left(\frac{m_{\text{jj}b}-m_{\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}}{0.1m_{\text{jj}b}}\right)^{2}}\,\right]\,,$

where $m_{W}=$80.4\text{\,}\mathrm{GeV}$$ and
$m_{t}=$172.5\text{\,}\mathrm{GeV}$$ are the nominal
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
boson and top quark masses, and $m_{\text{jj}}$ and $m_{\text{jj}b}$ are the
invariant masses of
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
boson and top quark candidates formed from jet combinations in each event. The
“minimum” refers to the minimum value from all possible jet combinations (of
one
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jet and two additional untagged jets) that would give a
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
boson candidate and a corresponding top candidate. The factor of $0.1$ in the
denominators is chosen to approximate the experimental dijet mass resolution.
The
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
boson candidates are formed from any pair of central jets in the event and the
top quark candidates are then reconstructed by pairing the
$\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
boson candidates with any remaining
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
Higgs boson candidate jets. The
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
discriminant is designed to quantify the likelihood that an event contains a
hadronic top quark decay. Events with
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<1.5$
are rejected. This reduces the
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
background by a factor of about 2 in simulated events, for a small loss of
signal efficiency, of around 15%, and a similar reduction in the
non-$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$,
multijet background.

In order to further reduce the overall background contamination, events in the
ggF signal region are also required to have reconstructed Higgs bosons that
satisfy a pseudorapidity separation
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$
. No such requirement is imposed in the VBF signal region, since SM VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal events tend to have a larger
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$.

A final analysis selection criterion to test the compatibility of events with
the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
decay is applied in both the ggF and VBF selections. A discriminant
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
is defined as:

$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=\sqrt{\left(\frac{m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}-$124\text{\,}\mathrm{GeV}$}{0.1\
m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}}\right)^{2}+\left(\frac{m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}-$117\text{\,}\mathrm{GeV}$}{0.1\
m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}}\right)^{2}}\,,$

where $m_{H1}$ and $m_{H2}$ are the masses of the leading and subleading
reconstructed Higgs boson candidates respectively. The values of
$124\text{\,}\mathrm{GeV}$ and $117\text{\,}\mathrm{GeV}$ in the
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
definition are chosen in accord with the centers of the $m_{H1}$ and $m_{H2}$
distributions for correctly paired signal events from simulation. Events are
required to have
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<1.6$
to be included in the Signal Region (SR) of the analysis.

Both the ggF and VBF signal regions are subdivided into a number of orthogonal
categories in order to better isolate the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal and improve the analysis sensitivity. The
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
and
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$
quantities are used to define six orthogonal ggF categories. The categories
are defined by two intervals in
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
with boundaries at 0, 0.95, and 1.6, and three in
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$,
with boundaries at 0, 0.5, 1, and 1.5 . In the VBF signal region, two
categories are defined using the
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$
quantity, with the dividing boundary at
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|=1.5$
. The
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$
category is more sensitive to VBF signals with non-SM couplings, while the
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.5$
category is more sensitive to SM VBF production.

The reconstructed invariant mass of the Higgs boson candidate pair,
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
is used as the discriminating variable for all analysis regions and categories
when extracting results, as detailed in Section 8. The
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution is found to have significant separation power between background
and signal, for all the different values of coupling modifiers. The binning of
the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions may vary between categories and is chosen in order to both
maintain discrimination power and limit the expected statistical uncertainty
in each bin to less than approximately $30\%$. This 30% limit ensures that the
assumptions used in the statistical procedure, outlined in Section 8, are
satisfied. In the VBF signal region, only events with
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>$400\text{\,}\mathrm{GeV}$$
are considered, as the background in the lower
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
region was found to be inadequately modeled by the data-driven method
described in Section 6 when applied to the
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
control data sample (also described in Section 6). For the ggF signal region,
no requirements on
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
are applied.

All the selection steps of the analysis are summarized in Figure 3. The yields
in the data and the simulated signal samples for some typical coupling values
are shown in Table 3. This sample of data events is referred to as $4b$ events
hereafter.

Figure 3: A flowchart summarizing the nine selection criteria used for the VBF and ggF analysis selections. Events must satisfy selection criteria 1–3 in order to be considered for either analysis signal region. Events failing to satisfy any of the selection criteria 4–6 are considered for inclusion in the ggF signal region, while those satisfying selection criteria 4–6 are considered for the VBF signal region.  Table 3: The yields of data and various example ggF and VBF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ signal models at each step of the analysis selection. The “Preselection” entry denotes an initial selection requiring at least four jets with $p_{\text{T}}>$40\text{\,}\mathrm{GeV}$$, at least two of which are $b$-tagged. Events which satisfy the “VBF selection” requirements are considered as part of the VBF signal region of the analysis, while the rest are considered for the ggF signal region. The signal yields are taken from simulation and are normalized by their theoretical cross-sections and the integrated luminosity of $126\text{\,}{\mathrm{fb}}^{-1}$. Corrections for differences in the $b$-tagging efficiency and trigger acceptance between data and simulation are applied starting from the “Trigger class” requirement.  |  Data | ggF Signal | VBF Signal
---|---|---|---
|  | SM | $\kappa_{\lambda}=10$ | SM | $\kappa_{2V}=0$
Common preselection |
Preselection | $5.70\text{\times}{10}^{8}$ | $530$ | $7300$ | $22$ | $630$
Trigger class | $2.49\text{\times}{10}^{8}$ | $380$ | $5300$ | $16$ | $410$
ggF selection |
Fail VBF selection | $2.46\text{\times}{10}^{8}$ | $380$ | $5200$ | $14$ | $330$
At least 4 $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged central jets | $1.89\text{\times}{10}^{6}$ | $86$ | $1000$ | $1.9$ | $65$
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$ | $1.03\text{\times}{10}^{6}$ | $72$ | $850$ | $0.94$ | $46$
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>1.5$ | $7.51\text{\times}{10}^{5}$ | $60$ | $570$ | $0.74$ | $43$
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<1.6$ (ggF signal region) | $1.62\text{\times}{10}^{4}$ | $29$ | $180$ | $0.24$ | $23$
VBF selection |
Pass VBF selection | $3.30\text{\times}{10}^{6}$ | $5.2$ | $81$ | $2.2$ | $71$
At least 4 $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged central jets | $2.71\text{\times}{10}^{4}$ | $1.1$ | $15$ | $0.74$ | $28$
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>1.5$ | $2.18\text{\times}{10}^{4}$ | $1.0$ | $11$ | $0.67$ | $26$
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<1.6$ | $5.02\text{\times}{10}^{2}$ | $0.48$ | $3.1$ | $0.33$ | $17$
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>$ $400\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ (VBF signal region) | $3.57\text{\times}{10}^{2}$ | $0.43$ | $1.8$ | $0.30$ | $16$

##  6 Background modeling

After the selection described above, about 90% of the background events come
from multijet processes (excluding top quark production), with the
approximately 10% remainder almost entirely composed of
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
events. This background composition was determined by applying the full event
selection to simulated samples of the various processes and comparing the
yields with the total background estimate in the SR; it is purely meant to be
indicative and is not used for deriving any results. The background is modeled
using the fully data-driven technique described below.

The background estimation makes use of an alternative set of events, which
pass the same
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-jet
triggers and satisfy all the same selection criteria as the
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events, with one difference: they are required to contain exactly two
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets. This sample, referred to hereafter as
“$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$”,
has about two orders of magnitude more events than the
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
sample, hence the presence of any
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
signal in it is negligible, making it suitable for the background estimation.
The jets selected to form the two Higgs boson candidates in the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events are the two
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets and the two untagged central jets with the highest $p_{\text{T}}$
(excluding the VBF jets in the VBF categories).

The kinematic properties of the $2b$ and $4b$ events are not expected to be
identical, partly due to different processes contributing to the two samples,
but also due to differences in the trigger acceptance and because the
probability of tagging a
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-jet
varies as a function of jet $p_{\text{T}}$ and $\eta$. Therefore, a
reweighting function is required, which, when applied to the $2b$ events, maps
their kinematic distributions onto the corresponding $4b$ distributions. This
function is derived using the $2b$ and $4b$ events in a Control Region (CR)
surrounding the SR in the reconstructed ($m_{H1}$, $m_{H2}$) plane and then
applied to the $2b$ events in the SR to produce the background estimate. The
“inner edge” of the CR is defined by
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=1.6$
and the “outer edge” by the circle:

$R_{\text{CR}}=\sqrt{\left(m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}-1.05\cdot$124\text{\,}\mathrm{GeV}$\right)^{2}+\left(m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}-1.05\cdot$117\text{\,}\mathrm{GeV}$\right)^{2}}=$45\text{\,}\mathrm{GeV}$\,.$

The shift of the center of the above circle by a factor of 1.05, relative to
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}=0$,
is found to be the optimal trade-off between having a good number of events
outside of the SR and avoiding the low $m_{H1}$/$m_{H2}$ regions, where the
differences between $2b$ and $4b$ kinematic distributions are larger. The CR
is split into four roughly equal directional quadrants, defined by
$45^{\circ}$ and $135^{\circ}$ lines passing through the SR center, (124, 117)
GeV. The four quadrants are given labels based on compass directions: the
upper quadrant $\text{Q}_{\text{N}}$, the lower $\text{Q}_{\text{S}}$, the
left $\text{Q}_{\text{W}}$, and the right $\text{Q}_{\text{E}}$. The above
lines also define four quadrants, with the same names as above, in the SR.
Events in CR $\text{Q}_{\text{N}}$ and $\text{Q}_{\text{S}}$, hereafter
referred to as CR1, are used to derive the reweighting function for the
nominal background estimate, while an alternative reweighting function,
derived from the CR events in $\text{Q}_{\text{E}}$ and $\text{Q}_{\text{W}}$
(referred to hereafter as CR2) is used to define a systematic uncertainty
related to the reweighting function interpolation into the SR, as detailed in
Section 7. The boundaries of the SR, CR1, and CR2 in the reconstructed
($m_{H1}$, $m_{H2}$) plane are shown in Figure 4. The horizontal and vertical
bands of lower event density around $80\text{\,}\mathrm{GeV}$ visible in these
plots are caused by the
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
selection criterion.

((a))  ((b))

Figure 4: The mass planes of the reconstructed Higgs boson candidates for the
(a) ggF and (b) VBF signal regions of the analysis, shown for the
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data events. In (a), the analysis selection up to step 8 (as outlined in
Figure 3) of the ggF selection has been applied, while in (b), the analysis
selection up to step 7 of the VBF selection has been applied. The continuous
red line describes the Signal Region (SR), the dashed line describes Control
Region 1 (CR1) and the dotted line describes Control Region 2 (CR2).

The reweighting function has the form:

$w(\vec{x})=\frac{p_{4b}(\vec{x})}{p_{2b}(\vec{x})},$ (1)

where $p_{4b}(\vec{x})$ and $p_{2b}(\vec{x})$ are the probability density
functions for $4b$ and $2b$ data, respectively, over a set of kinematic
variables $\vec{x}$. The computation of $w(\vec{x})$ is a density ratio
estimation problem, for which a variety of approaches exist. The method
employed in this analysis is modified from Refs. [90, 91] and makes use of an
artificial neural network (NN). This NN is trained on $2b$ and $4b$ CR1 data
(or CR2 data, for determining systematic uncertainties, as described Section
7). The training minimizes the following loss function:

$\mathcal{L}(w(\vec{x}))=\int
d\vec{x}\left[\sqrt{w(\vec{x})}p_{2b}(\vec{x})+\frac{1}{\sqrt{w(\vec{x})}}p_{4b}(\vec{x})\right]\,,$

The function in Eq. (1) minimizes this loss by equalizing the contributions
from the two terms. The kinematic variables used to make up $\vec{x}$ are
listed in Table 4 for the ggF and VBF signal regions; they are among those
kinematic variables that exhibit larger differences between the $2b$ and $4b$
events. The NN used in the ggF signal region has three densely connected
hidden layers of 50 nodes, each with a rectified linear unit activation
function [92], and a single-node linear output. A similar architecture is
chosen for the NN used in the VBF signal region, except that only 20 nodes are
used in each of the three hidden layers. This reflects the fact that the $2b$
and $4b$ sample sizes in the VBF signal region are nearly two orders of
magnitude smaller than the corresponding ones in the ggF signal region. This
is also the reason behind the choice to perform the NN training in the VBF
signal region for all data-taking years together, with the year index as a
one-hot encoded input feature.555One-hot encoding is a standard technique in
machine learning. For example, for the data-taking years in the VBF
reweighting, instead of presenting the year numbers as input features to the
NN, one-hot encoding uses three input features: $(1,0,0)$ for 2016, $(0,1,0)$
for 2017, and $(0,0,1)$ for 2018. For the ggF signal region, a dedicated
reweighting is derived for each year separately, which, thanks to the adequate
sample sizes, deals better with the different levels of disparity between
2$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
and
4$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
distributions, due to the differences in the trigger conditions from year to
year. Finally, in order to ensure that there are adequate numbers of
4$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events for both the ggF and VBF NN trainings, these trainings are performed
inclusively, before separating the events into the
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$
categories (the CR events have
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>1.6$,
hence it would not be possible to separate them into the
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
categories defined for the SR events). Both
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$
and
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
are found to be insensitive to the kinematic reweighting, and so the inclusive
training is not expected to introduce any additional bias when separating the
events into the various categories.

Table 4:  The set of input variables used for the 2$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ to 4$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ reweighting in the ggF and VBF channels respectively.  ggF | VBF
---|---
1. $\log(p_{\mathrm{T}})$ of the $2^{\mathrm{nd}}$ leading Higgs boson candidate jet 2. $\log(p_{\mathrm{T}})$ of the $4^{\mathrm{th}}$ leading Higgs boson candidate jet 3. $\log(\Delta R)$ between the closest two Higgs boson candidate jets 4. $\log(\Delta R)$ between the other two Higgs boson candidate jets 5. Average absolute $\eta$ value of the Higgs boson candidate jets 6. $\log(p_{\mathrm{T}}$) of the di-Higgs system 7. $\Delta R$ between the two Higgs boson candidates 8. $\Delta\phi$ between jets in the leading Higgs boson candidate 9. $\Delta\phi$ between jets in the subleading Higgs boson candidate 10. $\log(X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}})$ 11. Number of jets in the event 12. Trigger class index as one-hot encoder |  1. Maximum dijet mass from the possible pairings of the four Higgs boson candidate jets 2. Minimum dijet mass from the possible pairings of the four Higgs boson candidate jets 3. Energy of the leading Higgs boson candidate 4. Energy of the subleading Higgs boson candidate 5. Second-smallest $\Delta R$ between the jets in the leading Higgs boson candidate (from the three possible pairings for the leading Higgs candidate) 6. Average absolute $\eta$ value of the four Higgs boson candidate jets 7. $\log(X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}})$ 8. Trigger class index as one-hot encoder 9. Year index as one-hot encoder (for years inclusive training)

((a))  ((b))

Figure 5:  Comparison of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
(yellow histogram with hatching) and
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
(black points with error bars)
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions, for events in Control Region 1 (CR1) of the ggF signal region
from the 2018 data: (a) before the kinematic reweighting of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events, with only a normalization factor applied; and (b) after the kinematic
reweighting of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events. The error bars indicate the statistical uncertainty of the
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, while the hatching indicates the statistical uncertainty of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data. The latter is only the Poisson uncertainty of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, in (a), while in (b), it also includes the uncertainty from the
bootstrap procedure described in Section 7. The hatching in (a) is narrower
than the line width of the plotted histogram.

((a))  ((b))

Figure 6:  Comparison of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
(yellow histogram with hatching) and
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
(black points with error bars)
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions, for events in Control Region 1 (CR1) of the VBF signal region:
(a) before the kinematic reweighting of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events, with only a normalization factor applied; and (b) after the kinematic
reweighting of the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
events. The error bars indicate the statistical uncertainty on the
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, while the hatching indicates the statistical uncertainty on the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data. The latter is only the Poisson uncertainty on the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, in (a), while in (b), it also includes the uncertainty from the
bootstrap procedure described in Section 7.

In order to estimate and mitigate the impact of the varying initial conditions
and limited size of the training samples on the NN training, the deep
ensembles technique [93] is used together with a bootstrap resampling [94] of
the training data. This entails constructing a set of training datasets by
sampling with replacement from the original dataset. In this analysis, this is
approximated by the usage of different random training weights, following a
Poisson distribution with $\mu=1$, for each event in each training. The NN is
trained independently on each element of this set, using different initial
conditions each time. This results in an ensemble of reweighting functions.
Each reweighting function is further multiplied by a normalization factor,
such that the number of reweighted $2b$ events is equal to the number of $4b$
events in the region where the NN is trained. In this analysis, the ensembles
comprise 100 reweighting functions each, hence 100 weights are calculated for
each $2b$ event in the SR. The background estimate uses the mean of these
weights for each event, and the variation of the background predictions from
the ensemble of reweighting functions is used to estimate a systematic
uncertainty for the stability of the NN training procedure, as described in
Section 7.

The effect of the above reweighting procedure in CR1, where the reweighting
function is derived, is illustrated in Figure 5 for the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution of the ggF-selected events and in Figure 6 for the
$X_{\mathit{{W}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution of the VBF-selected events. The reweighted
“$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$”
distributions agree with the corresponding
“$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$”
distributions to within about 10% for most of the phase space, with some
larger deviations observed in bins near the tails of the distributions where
fewer data events are available. A large number of additional kinematic
variables were also studied before and after applying the reweighting in order
to validate the performance of the NN. For all variables, the level of
agreement, as quantified by the $\chi^{2}$ metric, either improves after the
reweighting or, for variables where the
“$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$”
and
“$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$”
distributions are already similar, changes only slightly.

The background modeling procedure was tested and found to produce good results
in a large simulated
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
sample and a much smaller sample of simulated
(non-$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$)
multijet events in the SR. The procedure was also tested in several control
data samples orthogonal to the nominal event selection, where the presence of
any
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal is negligible and the $4b$ events in the corresponding SR can be
compared with the reweighted SR $2b$ events without any bias. These samples,
summarized in Table 5, include: (a) events satisfying all the $2b$/$4b$ ggF
selection criteria, with the difference that the
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$
cut is inverted; (b) events satisfying all the $2b$/$4b$ selection criteria,
except that the center of the SR (and hence also of CR1 and CR2) is shifted,
to avoid any overlap with the nominal SR; and (c) events that satisfy all the
same $4b$ selection criteria, except that, in terms of
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagging,
they contain exactly three
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged
jets, and the fourth jet is taken as the highest-$p_{\text{T}}$ jet that fails
a looser working point of the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagging
algorithm (one that gives 85% efficiency for
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-jets
in simulated
$\mathit{{t}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{t}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
events). This sample, hereafter referred to as
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$,
has about one order of magnitude more events than the $4b$ sample and a
negligible amount of
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal; hence it is used to derive a nonclosure systematic uncertainty for the
reweighting procedure, as discussed in Section 7. No significant background
modeling nonclosure was observed in the other control data samples.

Table 5:  A summary of all the data samples used in the analysis. For each control sample, the variations from the nominal analysis are noted explicitly.  Data Sample | Definition | Usage
---|---|---
Signal Region (SR) |  Events with $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<1.6$ |  Defines signal region in the $m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}$–$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}$ plane
Control Region (CR) |  Events with $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>1.6$ and $R_{\text{CR}}<45$ GeV |  Defines control region in the $m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}$–$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}$ plane for background estimation (ggF and VBF)
Shifted validation regions |  Shift the center of the SR in the $m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1}$–$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}2}$ plane to avoid overlap with the nominal SR |  Background estimation validation (ggF only)
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ |  Events with at least 4 $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged central jets |  Final analysis sample (ggF and VBF)
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ |  Events with exactly 2 $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged central jets plus at least two additional untagged central jets |  Background estimation (ggF and VBF)
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$ |  Events with exactly 3 $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagged central jets plus at least one central jet failing a looser $\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$-tagging requirement |  Background estimation validation (ggF and VBF), additional background modeling uncertainty (ggF only)
Reverse $|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$ |  $2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ and $4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ events with $|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.5$ |  Background estimation validation (ggF only)

##  7 Systematic uncertainties

The uncertainties with the greatest impact on the analysis sensitivity are
those arising from the data-driven background estimate described in Section 6.
These uncertainties have two main sources: the limited sample sizes in the CR
and SR, and physical differences between the CR, where the
2$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
reweighting function is derived, and the SR, where it is applied.

As described in Section 6, the ensemble of 100 reweighting functions results
in 100 separate background predictions. An
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
histogram can be constructed from each of these predictions, and the standard
deviation of the predictions in each bin is taken as the bootstrap
uncertainty. The uncertainty is treated as uncorrelated across
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
bins.

An additional statistical uncertainty results from the limited sample size of
the
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
SR dataset in which the trained background reweighting network is applied to
obtain the final background estimate. A Poisson uncertainty is taken for each
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
bin, which is combined in quadrature with the bootstrap uncertainty described
above.

For the background estimate, the uncertainty component related to the
kinematic differences between the SR and CR1 is evaluated by using alternative
predictions from the CR2 region. Four alternative background estimates are
produced by applying the CR1-derived weights to three of the SR quadrants, and
CR2-derived weights to the one remaining SR quadrant. For example, one
alternative background estimate is obtained by applying CR1-derived weights to
Q${}_{\text{S}}$, Q${}_{\text{E}}$ and Q${}_{\text{W}}$, and CR2-derived
weights to Q${}_{\text{N}}$. Each of these four background predictions is
symmetrized around the nominal
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution to construct a two-sided uncertainty. Since the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution differs across the four SR quadrants, substituting the CR2-based
prediction for the CR1-based prediction in each of the four SR quadrants
separately and utilizing a four-component uncertainty gives the fit model
greater flexibility to describe these
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
variations with finer granularity. In the ggF signal region, these
uncertainties are taken to be uncorrelated across the datasets from the three
different years. In both the ggF and VBF signal regions, the uncertainty is
treated as correlated across the analysis categories.

An additional closure uncertainty is estimated by applying the full background
modeling procedure to the
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
sample instead of the
4$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
sample. The predicted
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distribution in the various analysis categories is then compared with the
observed
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
data in the SR. For the VBF signal region, no statistically significant
difference between the prediction and observation is found, and hence no
additional uncertainty is applied. For the ggF signal region, an additional
uncertainty is evaluated in each category from the observed differences
between the predicted and observed
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions. For
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
bins in which the predicted and observed values differ by less than $1\sigma$,
where $\sigma$ is obtained from all other background modeling uncertainties
combined, no additional uncertainty is applied. For
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
bins where the predicted and observed values differ by more than $1\sigma$,
the amount beyond $1\sigma$ is averaged with the corresponding amounts in the
two adjacent bins, to limit the impact of statistical fluctuations, and is
symmetrized around the nominal prediction to construct a two-sided
uncertainty. This nonclosure uncertainty has a much smaller impact on the
analysis sensitivity than the other sources of background modeling
uncertainty.

Several detector modeling uncertainties are evaluated and included. These
affect only the signal description, as the background is estimated entirely
from data. Uncertainties in the jet energy scale and resolution, as well as
the JVT, are treated according to the prescription in Refs. [80, 81].
Additional uncertainties arising from the correction of the simulated pileup
distribution are treated according to the prescription in Ref. [95].
Uncertainties in the $b$-tagging efficiency are treated according to the
prescription in Ref. [96]. Uncertainties in the trigger efficiencies are
evaluated from measurements of per-jet online efficiencies for both jet
reconstruction and $b$-tagging, which are used to compute event-level
uncertainties. These are then applied to the simulated events as overall
weight variations. The uncertainty in the integrated luminosity used in this
analysis is in the range 2.0%–2.4% for the three years of data taking and 1.7%
for the entire dataset [97], obtained using the LUCID-2 detector for the
primary luminosity measurements [98].

Several sources of theoretical uncertainty affecting the signal models are
considered as described below. Uncertainties due to modeling of the parton
shower and underlying event are evaluated by comparing results between two
generators for these parts of the calculation: the nominal Pythia 8 and the
alternative Herwig 7\. This is found to have an effect of roughly $10\%$ on
the ggF and VBF signal acceptances, and a negligible impact on the shape of
the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions. The parton showering uncertainty is derived within each
analysis SR category; the uncertainty is observed to reach approximately
$40\%$ for a given production mode in some categories in which the acceptance
is small for that mode. Uncertainties in the matrix element calculation are
evaluated by raising and lowering the factorization and renormalization scales
used in the generator by a factor of two, both independently and
simultaneously. This results in an effect of typically $2\%$ for both ggF and
VBF, with a maximum effect of about $6\%$ in certain analysis categories. PDF
uncertainties are evaluated using the PDF4LHC_NLO_MC set [54] by calculating
the signal acceptance for each replica and taking the standard deviation. The
magnitude of this uncertainty is typically found to be less than $1\%$ in both
the ggF and VBF signal acceptances, with a maximum magnitude of approximately
$2\%$. Theoretical uncertainties in the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}$
branching ratio [14] are included, amounting to an approximately $3.5\%$
overall uncertainty in the signal normalization. The dependence of the
branching ratio uncertainty on $\kappa_{\lambda}$ is neglected. Theoretical
uncertainties in the ggF and VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
cross-sections arising from uncertainties in the PDF and $\alpha_{\text{s}}$,
as well as the choice of renormalization scheme and the scale of the top
quark’s mass, are taken from Refs. [14, 99, 11]. The cross-section
uncertainties are included in the derivation of the upper limits on the ggF,
VBF, and combined
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal strengths, as well as the likelihood-based constraints on the values of
the $\kappa_{\lambda}$ and $\kappa_{2V}$ modifiers, as presented in Section 8.

An additional signal modeling systematic uncertainty is evaluated for the
SMEFT and HEFT measurements. The
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
spectra of reweighted SMEFT/HEFT signal samples are compared against
explicitly generated samples for a select number of coefficient variations. A
two-component normalization uncertainty is derived by taking the average of
the relative deviations across the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
bins in the ranges of
$280<m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<$936\text{\,}\mathrm{GeV}$$
and
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>$936\text{\,}\mathrm{GeV}$$.
The use of separate components in the low- and
high-$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
regions prevents the level of agreement in the more populated
low-$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
region from overconstraining the uncertainty in the more sparsely populated
high-$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
region.

##  8 Results

The analysis results are obtained using a maximum-likelihood fit performed in
bins of reconstructed
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$.
For the ggF signal region, the fit is performed simultaneously across the
different data-taking years (2016, 2017, and 2018), while for the VBF signal
region, the fit is performed inclusively on the data from all years.

The likelihood function used to construct the test statistic has a standard
form, consisting of a product of Poisson distributions for the yields in each
bin and constraint functions for nuisance parameters describing systematic
uncertainties. For uncertainties due to the limited sample size in data or
simulation, the constraint is a Poisson distribution. For all other systematic
uncertainties, the constraint is a Gaussian distribution. Where systematic
uncertainties are deemed to be uncorrelated, independent nuisance parameters
are introduced. Uncertainties in the luminosity and signal modeling are
treated as fully correlated between the analysis categories and, for ggF, the
data-taking years. Each component of the quadrant-derived uncertainty covering
the kinematic differences between the SR and CR1 regions is correlated across
the data-taking years for the ggF region. The components are correlated across
analysis categories within the ggF and VBF signal regions, but not between the
ggF and VBF signal regions. All other uncertainties in the background model
are treated as uncorrelated across the different categories and data-taking
years. The statistical model is implemented using RooFit [100].

The hypothesis of the presence of a signal is tested using the profile
likelihood ratio [101]. The signal strength of the combined ggF and VBF signal
process, $\mu_{\text{ggF}+\text{VBF}}$
($=\sigma_{\text{ggF}+\text{VBF}}/\sigma_{\text{ggF}+\text{VBF}}^{\text{SM}}$),
is chosen as the parameter of interest (POI) and is a free parameter in the
fit. The relative contributions of the ggF and VBF signals to the total signal
model are fixed to their predicted values. The profile likelihood ratio takes
the following form:

$-2\Delta\ln\lambda(\mu)=-2\ln\Bigg{(}\frac{L(\mu,\hat{\hat{\bm{\theta}}})}{L(\hat{\mu},\hat{\bm{\theta}})}\Bigg{)}\,,$

where $\mu$ is the POI and $\bm{\theta}$ represents the nuisance parameters.
The numerator represents the conditional maximum-likelihood fit, in which the
nuisance parameters are set to their profiled values $\hat{\hat{\bm{\theta}}}$
for which the likelihood is maximized for a fixed value of $\mu$. The
denominator represents the unconditional likelihood fit, where both $\mu$ and
$\bm{\theta}$ are set to the values which jointly maximize the likelihood,
$\hat{\mu}$ and $\hat{\bm{\theta}}$, respectively.

The observed distributions in
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$,
as well as the predicted background and example signal shapes, are presented
in Figure 7 for each of the six ggF categories (with all data-taking years
combined, for presentation purposes). The distributions of the expected
background are obtained using the best-fit values of the nuisance parameters
in the fit to the data with the background-only hypothesis. The corresponding
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
distributions in the two VBF categories are shown in Figure 8. The signal
shape for $\kappa_{2V}=0$ in Figure 8(a) clearly shows the impact of the
divergences in Figures 2(b) and 2(c) not canceling out as discussed in Section
1. While the deviations from the SM studied in this analysis are below the
level that violates unitarity, this behavior makes the VBF topology in this
analysis particularly sensitive to $\kappa_{2V}$. The observed number of data
events, predicted number of background events, and expected number of signal
events for the SM ggF and VBF signals are summarized for each of the analysis
categories in Table 6.

Table 6: The yields in each analysis category of the data, expected background, and expected SM ggF and VBF signals. The expected background yields are obtained using a fit to the data with the background-only hypothesis; the quoted uncertainties are the sum in quadrature of all the per-bin systematic uncertainties. The expected signal yields are obtained from simulation. Category | Data | Expected | ggF Signal | VBF Signal
---|---|---|---|---
|  | Background | SM | SM
ggF signal region |
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<0.5$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$ | $1940$ | $1935\pm 25$ | $7.0$ | $0.038$
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<0.5$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$ | $3602$ | $3618\pm 37$ | $6.5$ | $0.036$
$0.5<|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.0$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$ | $1924$ | $1874\pm 21$ | $5.1$ | $0.037$
$0.5<|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.0$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$ | $3540$ | $3492\pm 35$ | $4.7$ | $0.040$
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.0$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$ | $1880$ | $1739\pm 22$ | $2.9$ | $0.043$
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.0$, $X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$ | $3285$ | $3212\pm 37$ | $2.8$ | $0.041$
VBF signal region |
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$ | $116$ | $125.3\pm 4.4$ | $0.37$ | $0.090$
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.5$ | $241$ | $230.6\pm 5.3$ | $0.06$ | $0.21$

((a))  ((b))  ((c))

((d))  ((e))  ((f))

Figure 7:  Distributions of the reconstructed
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
in data (shown by the black points) and the estimated background (shown by the
yellow histograms), in each of the six
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
categories in the ggF signal region: (a)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<0.5$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$;
(b)
$0.5<|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.0$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$;
(c)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.0$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}<0.95$;
(d)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<0.5$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$;
(e)
$0.5<|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.0$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$;
and (f)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.0$,
$X_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}>0.95$.
The contributions from the different data-taking years are combined in each
category for presentation purposes. The hatching shows the total uncertainty
of the background estimate. The distribution of the expected background is
obtained using the best-fit values of the nuisance parameters in the fit to
the data with the background-only hypothesis. Distributions of the SM and
$\kappa_{\lambda}=6$ signal models are overlaid, scaled so as to be visible on
the plot, and the scaling for each signal model is the same across the six
categories. The lower panels show the ratio of the observed data yield to the
predicted background in each bin. Events in the underflow and overflow bins
are counted in the yields of the initial and final bins respectively.

((a))  ((b))

Figure 8: Distributions of the reconstructed
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
in data (shown by the black points), the estimated background (shown by the
yellow histograms), in each of the two
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|$
categories in the VBF signal region: (a)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|<1.5$
and (b)
$|\Delta\eta_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}|>1.5$.
The hatching shows the total uncertainty of the background estimate. The
distribution of the expected background is obtained using the best-fit values
of the nuisance parameters in the fit to the data with the background-only
hypothesis. Distributions for three choices of couplings are shown: the SM,
$\kappa_{\lambda}=6$, and $\kappa_{2V}=0$ (with all other couplings set to
their SM values in the last two models), scaled so as to be visible on the
plot. The lower panels show the ratio of the observed data yield to the
predicted background in each bin. Events in the overflow bins are counted in
the yields of the final bins.

An upper limit on the combined ggF and VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal strength $\mu_{\text{ggF+VBF}}$ is computed using the asymptotic
formula [101] and based on the CL${}_{\text{s}}$ method [102]. The observed
(expected) 95% CL upper limit on $\mu_{\text{ggF+VBF}}$ is found to be 5.4
(8.1). The expected upper limits are obtained using a background-only
hypothesis, excluding a $HH$ signal. The upper limit on the combined
$\mu_{\text{ggF+VBF}}$, as well as upper limits on the individual
$\mu_{\text{ggF}}$ ($=\sigma_{\text{ggF}}/\sigma_{\text{ggF}}^{\text{SM}}$)
and $\mu_{\text{VBF}}$
($=\sigma_{\text{VBF}}/\sigma_{\text{VBF}}^{\text{SM}}$), are summarized in
Table 7. For the individual $\mu_{\text{ggF}}$ and $\mu_{\text{VBF}}$ limits,
the results are derived by treating the other production mode (VBF when
placing limits on $\mu_{\text{ggF}}$, and vice-versa) as a background process,
with its normalization only loosely constrained in the fit. Compared to the
previous ATLAS measurement of ggF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production in the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
decay channel (using $27\text{\,}{\mathrm{fb}}^{-1}$ of early Run 2 data)
[26], the upper limit on the ggF cross-section is over 50% lower, with
approximately 20% of this improvement arising from advances in analysis
techniques and object reconstruction. Similarly, compared to the previous
ATLAS measurement of VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
production in the
$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$$\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}$
decay channel, which used $126\text{\,}{\mathrm{fb}}^{-1}$ of data collected
between 2016 and 2018 [27], the upper limit on the VBF
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
cross-section is over 75% lower, with this improvement arising entirely from
advances in analysis technique and object reconstruction.

Table 7: The observed and expected upper limits on the SM ggF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ production cross-section $\sigma_{\text{ggF}}$, SM VBF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ production cross-section $\sigma_{\text{VBF}}$, and combined SM ggF and VBF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ production cross-section $\sigma_{\text{ggF+VBF}}$ at the 95% CL, expressed as multiples of the corresponding SM cross-sections. The expected values are shown with corresponding one- and two-standard-deviation error bounds, and they are obtained using a background-only fit to the data. When extracting the limits on $\sigma_{\text{ggF+VBF}}$, the relative contributions of ggF and VBF production to the total cross-section are fixed to the SM prediction. | Observed Limit | $-2\sigma$ | $-1\sigma$ | Expected Limit | $+1\sigma$ | $+2\sigma$
---|---|---|---|---|---|---
$\mu_{\text{ggF}}$ | $5.5$ | $4.4$ | $5.9$ | $8.2$ | $12.4$ | $19.6$
$\mu_{\text{VBF}}$ | $130$ | $70$ | $100$ | $130$ | $190$ | $280$
$\mu_{\text{ggF}+\text{VBF}}$ | $5.4$ | $4.3$ | $5.8$ | $8.1$ | $12.2$ | $19.1$

The total uncertainty in the upper limit of the cross-section is dominated by
the uncertainty sources related to the background modeling procedure and
theoretical predictions. With only the statistical uncertainties of the
reweighted
$2\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, observed
$4\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
data, and simulated signal samples included in the fit, the expected upper
limit on $\mu_{\text{ggF+VBF}}$ is found to be $6.0$ times the SM prediction.
Including the uncertainty sources resulting from the background estimation
(the bootstrap uncertainty, the uncertainty from the kinematic differences
between the SR and CR1, and, in the ggF signal region, the
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$
nonclosure uncertainty), the expected upper limit on $\mu_{\text{ggF+VBF}}$ is
relaxed to $7.1$ times the SM prediction. The further reduction of sensitivity
to the value of 8.1, as quoted in Table 7, is driven primarily by the
uncertainties arising from theoretical predictions. The relative impact of the
various sources of systematic uncertainty on the expected upper limit on
$\mu_{\text{ggF+VBF}}$ is summarized in Table 8.

Table 8: Breakdown of the dominant systematic uncertainties. The impact of the uncertainties on the expected upper limit on $\mu_{\text{ggF+VBF}}$ when re-evaluating the profile likelihood ratio after fixing the nuisance parameter(s) in question to its (their) best-fit value(s), while all remaining nuisance parameters remain free to float. The impact is shown in %. Only (groups of) systematic uncertainties that have an impact of at least 1% are shown. The impact of each experimental source of systematic uncertainty described in the text, as well as of all of them together, is less than 1%.  Source of Uncertainty | $\Delta\mu/\mu$
---|---
Theory uncertainties |
Theory uncertainty in signal cross-section |  $-9.0$%
All other theory uncertainties |  $-1.4$%
Background modeling uncertainties |
Bootstrap uncertainty |  $-7.1$%
CR to SR extrapolation uncertainty |  $-7.5$%
$3\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}1\text{f}$ nonclosure uncertainty |  $-2.0$%

Constraints are placed on the $\kappa_{\lambda}$ and $\kappa_{2V}$ modifiers
using two different interpretations, the first named the “$95\%$ CL” method
and the second named the “profile likelihood ratio” method. The former uses
the signal strength $\mu$ as the POI, while the latter uses the vector of
coupling modifiers $\bm{\kappa}=(\kappa_{\lambda},\kappa_{2V})$. The $95\%$ CL
method allows for interpretation as a traditional search for an arbitrarily
normalized set of signals with different shapes against an estimated
background, while the profile likelihood ratio method allows for
interpretation as to whether the data are compatible with the specific cross-
section and shape predictions of the $\kappa$ framework. The $95\%$ CL results
presented here offer a consistent comparison with previous ATLAS
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
measurements. The constraints obtained from the two interpretations are not
expected to be identical, as the two strategies employ slightly different
physical assumptions. In the profile likelihood ratio interpretation, the
signal strength is fixed to the prediction obtained for a specific coupling
modifier configuration, while for the $95\%$ CL interpretation, the signal
strength is allowed to float. The profile likelihood ratio method utilizes a
hypothesis consisting of the predicted background plus the SM
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal, while the $95\%$ CL results utilize a hypothesis containing only the
predicted background and no
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal. Given the relatively small size of the SM
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal compared to the predicted background, the use of different hypotheses
is not expected to have a significant effect. Additionally, $2\sigma$-level
constraints are quoted from the profile likelihood ratio interpretation, as
opposed to $95\%$ CL constraints.

The $95\%$ CL constraints on $\kappa_{\lambda}$ and $\kappa_{2V}$ are obtained
by determining the 95% CL upper limits on the cross-section as a function of
these coupling modifiers,
$\mu_{\text{ggF}+\text{VBF}}(\kappa_{\lambda},\kappa_{2V})$. Values of the
coupling modifiers $(\kappa_{\lambda},\kappa_{2V})$ are excluded if the
predicted cross-section of the signal model obtained with that configuration
is excluded at the 95% CL. The
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\rightarrow\mathit{{b}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}{}\mathit{{\mkern
0.0mu\overline{\mkern
0.0mu{b}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}{}$
branching ratio is fixed to the SM prediction in the likelihood fit and any
dependence on $\kappa_{\lambda}$ is ignored. Upper limits on the
$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
signal strength as a function of $\kappa_{\lambda}$ and $\kappa_{2V}$ are
shown in Figure 9, and the exclusion boundaries are summarized in Table 9.
With the values of the other modifiers ($\kappa_{V}$ and either $\kappa_{2V}$
or $\kappa_{\lambda}$, respectively) fixed to their SM value of 1, values of
$\kappa_{\lambda}$ beyond $[-3.9,11.1]$ and values of $\kappa_{2V}$ beyond
$[-0.03,2.11]$ are excluded.

((a))  ((b))

Figure 9: The observed 95% CL exclusion limits as a function of (a) $\kappa_{\lambda}$ (obtained using the signal strength $\mu_{\text{ggF}+\text{VBF}}$ as the POI) and (b) $\kappa_{2V}$ (obtained using the signal strength $\mu_{\text{VBF}}$ as the POI) from the combined ggF and VBF signal regions, as shown by the solid black line. In each case, the value of the other modifier is fixed to 1. The blue and yellow bands show respectively the $1\sigma$ and $2\sigma$ bands around the expected exclusion limits, which are shown by the dashed black line. The expected exclusion limits are obtained using a fit to the data with the background-only hypothesis. The dark red line shows in (a) the predicted combined ggF and VBF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ cross-section as a function of $\kappa_{\lambda}$ and in (b) the predicted VBF $\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$$\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$ cross-section as a function of $\kappa_{2V}$. The dark pink bands surrounding the predicted cross-section lines indicate the theoretical uncertainty of the cross-section, as taken from Ref. [99]. The band in (b) is smaller than the width of the plotted line.  Table 9: The observed and expected constraints on the $\kappa_{\lambda}$ and $\kappa_{2V}$ coupling modifiers at 95% CL. For each modifier, the constraints were extracted with all other modifiers fixed to the SM prediction. Parameter | Expected Constraint | Observed Constraint
---|---|---
| Lower | Upper | Lower | Upper
$\kappa_{\lambda}$ | $-4.6$ | $10.8$ | $-3.9$ | $11.1$
$\kappa_{2V}$ | $-0.05$ | $2.12$ | $-0.03$ | $2.11$

Figure 10 shows the $95\%$ CL exclusion limits in the two-dimensional plane of
the $\kappa_{\lambda}$–$\kappa_{2V}$ modifier space.

Figure 10: The observed $95\%$ CL exclusion limit obtained using the signal
strength $\mu_{\text{ggF}+\text{VBF}}$ as the POI in the two-dimensional
$\kappa_{\lambda}$ vs. $\kappa_{2V}$ space, obtained from the combined ggF and
VBF signal model, as shown by the solid black line. The blue and yellow bands
show respectively the $1\sigma$ and $2\sigma$ bands around the expected
exclusion limits, which are shown by the dashed black line. The star denotes
the SM prediction ($\kappa_{\lambda}=\kappa_{2V}=1$).

((a))  ((b))

Figure 11: The observed profile likelihood ratio scans for the (a)
$\kappa_{\lambda}$ and (b) $\kappa_{2V}$ coupling modifiers, shown by the
solid black line, using the coupling modifiers $\bm{\kappa}$ as the POIs. In
each case, the value of the other parameter is fixed to 1. The dashed blue
line shows the expected profile likelihood ratio, as obtained using a fit to
the data with the background-only hypothesis. The pink line indicates the
$2\sigma$ exclusion boundary.

((a))  ((b))

Figure 12:  (a) The observed profile likelihood ratio exclusion limits for the
two-dimensional $\kappa_{\lambda}$ vs. $\kappa_{2V}$ modifier space, shown by
the solid dark purple line at the $1\sigma$ level and the dashed turquoise
line at the $2\sigma$ level. The black cross denotes the best-fit values of
$(\kappa_{\lambda},\kappa_{2V})$. The expected exclusion limits are presented
in (b), where the solid pink line denotes the $1\sigma$-level exclusion and
the dashed orange line denotes the $2\sigma$-level exclusion. For both the
expected and observed limit plots, the black star indicates the SM prediction
($\kappa_{\lambda}=\kappa_{2V}=1$).

The alternative coupling modifier constraints are obtained using the profile
likelihood ratio interpretation, with the coupling modifiers
$\bm{\kappa}=(\kappa_{\lambda},\kappa_{2V})$ as the POIs, rather than the
signal strength $\mu$:

$-2\Delta\ln
L(\bm{\kappa})=-2\ln\Bigg{(}\frac{L(\bm{\kappa},\hat{\hat{\bm{\theta}}})}{L(\hat{\bm{\kappa}},\hat{\bm{\theta}})}\Bigg{)}\,.$

A scan of the profile likelihood ratio is performed as a function of the
coupling modifiers at discrete points to produce the curves shown in Figure
11. The best-fit value of $\kappa_{\lambda}$ is found to be 6.2 from the
profile likelihood scan. The observed pull of the best-fit $\kappa_{\lambda}$
value away from the SM value is due to a slight excess in the observed data in
the ggF signal region, specifically in the
low-$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$
range. The particular signal model in which $\kappa_{\lambda}$ is close to $6$
is favored due to a balance between two competing effects: the
$m_{\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{H}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}}$