HIG-21-002

HIG-21-002

# Search for Higgs boson pairs decaying to $\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$,
and $\Pgt\Pgt\Pgt\Pgt$ in proton-proton collisions at $\sqrt{s}=13\TeV$

###### Abstract

The results of a search for Higgs boson pair ($\PH\PH$) production in the
$\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$, and $\Pgt\Pgt\Pgt\Pgt$ decay modes are
presented. The search uses 138of proton-proton collision data recorded by the
CMS experiment at the LHC at a center-of-mass energy of 13from 2016 to 2018.
Analyzed events contain two, three, or four reconstructed leptons, including
electrons, muons, and hadronically decaying tau leptons. No evidence for a
signal is found in the data. Upper limits are set on the cross section for
nonresonant $\PH\PH$ production, as well as resonant production in which a new
heavy particle decays to a pair of Higgs bosons. For nonresonant production,
the observed (expected) upper limit on the cross section at 95% confidence
level () is 21.3 (19.4) times the standard model (SM) prediction. The observed
(expected) ratio of the trilinear Higgs boson self-coupling to its value in
the SM is constrained to be within the interval $-6.9$ to 11.1 ($-6.9$ to
11.7) at 95% , and limits are set on a variety of new-physics models using an
effective field theory approach. The observed (expected) limits on the cross
section for resonant $\PH\PH$ production amount to 0.18–0.90 (0.08–1.06) pb at
95% for new heavy-particle masses in the range 250–1000.

##  0.1 Introduction

Since the discovery of the Higgs () boson [1, 2, 3], many of its properties
have already been measured with high precision [4, 5, 6]. One important
property that remains largely unknown is the boson self-coupling. A precise
measurement of this coupling is necessary to determine the shape of the Higgs
potential, and thus verify that the mechanism breaking the electroweak gauge
symmetry is indeed the Higgs mechanism [7, 8, 9, 10, 11, 12] of the standard
model (SM) [13, 14, 15]. The SM predicts the existence of trilinear as well as
quartic boson self-couplings. Due to the very low cross section for triple
boson production, the quartic self-coupling will not be experimentally
accessible at the CERN LHC, even with the full integrated luminosity of
3000scheduled to be delivered after the high-luminosity LHC upgrade [16, 17].
The strength of the trilinear self-coupling, however, can be determined using
measurements of boson pair ($\PH\PH$) production.

In the SM, most $\PH\PH$ pairs are produced in two types of processes. The
Feynman diagrams for the dominant “gluon fusion” ($\Pg\Pg\PH\PH$) process at
leading order (LO) in perturbative quantum chromodynamics (QCD) are shown in
Fig. 1. The left “triangle” diagram amplitude varies proportionally to the
boson self-coupling ($\lambda$) and the Yukawa coupling of the top quark
($\text{y}_{\PQt}$), while the right “box” diagram amplitude is insensitive to
$\lambda$ and varies as $\text{y}_{\PQt}^{2}$. The triangle and box diagrams
interfere destructively, so the $\Pg\Pg\PH\PH$ cross section exhibits a strong
dependence on $\lambda$ and $\text{y}_{\PQt}$. The ratios of $\lambda$ and
$\text{y}_{\PQt}$ to their SM expectations are denoted as $\kappa_{\lambda}$
and $\kappa_{\PQt}$, respectively. By definition, these “coupling strength
modifiers” have values $\kappa_{\lambda}=1$ and $\kappa_{\PQt}=1$ in the SM.
The $\Pg\Pg\PH\PH$ cross section in the SM has been computed to be
$31.1^{+2.1}_{-7.2}$ fb at next-to-next-to-LO (NNLO) accuracy in QCD using the
FT${}_{\text{approx}}$ scheme, in which the true top quark mass is used for
the real radiation matrix elements, while the virtual part is computed using
an infinite top quark mass [18]. The SM cross section for the subdominant
“vector boson fusion” ($\PQq\PQq\PH\PH$) process has been computed at next-to-
NNLO accuracy in QCD and amounts to $1.73\pm 0.04$ fb [19].

Figure 1:  Leading order Feynman diagrams for SM nonresonant $\PH\PH$
production via gluon fusion, including the “triangle” diagram (left) and the
“box” diagram (right).

Deviations of these coupling strength modifiers from unity not only affect the
rate of $\PH\PH$ production, but also the kinematic distributions of the
$\PH\PH$ signal. The $\PH\PH$ invariant mass ($m_{\PH\PH}$) is particularly
sensitive to changes in $\kappa_{\lambda}$ and $\kappa_{\PQt}$, as these
couplings affect the triangle and box diagram amplitudes differently. Because
SM $\Pg\Pg\PH\PH$ and $\PQq\PQq\PH\PH$ production do not include a heavy
resonant particle, and typically result in a broad $m_{\PH\PH}$ distribution,
they are referred to as “nonresonant”. Changes in $\kappa_{\lambda}$ and
$\kappa_{\PQt}$ also influence the rate of single Higgs boson production as
well as the Higgs boson decay branching fractions [20, 21].

The presence of undiscovered particles or interactions, predicted by a variety
of theoretical models beyond the SM, may alter the $\PH\PH$ production rate as
well as observable kinematic distributions. Such particles could give rise to
loop diagrams similar to those shown in the upper row of Fig. 1. These
diagrams may significantly enhance the $\PH\PH$ production rate, as they occur
at the same loop level as $\PH\PH$ production in the SM. Since no particles
beyond those predicted by the SM have been observed so far, their mass may be
at the scale or higher, well above the scale of electroweak symmetry breaking.
Loop contributions of such heavy particles can be approximated as contact
interactions with the boson using an effective field theory (EFT) approach
[22, 23]. Following Ref. [24], the contact interactions relevant for $\PH\PH$
production are parametrized by the couplings $\text{c}_{\Pg}$,
$\text{c}_{2\Pg}$, and $\text{c}_{2}$, referring to the interactions between
two gluons and one boson, two gluons and two bosons, and two top quarks and
two bosons, respectively. The corresponding Feynman diagrams for
$\Pg\Pg\PH\PH$ production are shown in Fig. 2. The LO diagrams for
$\PQq\PQq\PH\PH$ production contain no gluons or top quarks, so the impacts of
$\text{c}_{\Pg}$, $\text{c}_{2\Pg}$, and $\text{c}_{2}$ are only considered in
the $\Pg\Pg\PH\PH$ signal in this publication.

Figure 2:  Leading order Feynman diagrams for nonresonant $\PH\PH$ production
via gluon fusion in an EFT approach, where loop-mediated contact interactions
between (left) two gluons and one boson, (middle) two gluons and two bosons,
and (right) two top quarks and two bosons are parametrized by three effective
couplings: $\text{c}_{\Pg}$, $\text{c}_{2\Pg}$, and $\text{c}_{2}$.

An excess of $\PH\PH$ signal events may also result from decays of new heavy
particles, denoted as X, into pairs of bosons. Various theoretical models of
new physics postulate such decays, in particular two-Higgs-doublet models [25,
26], composite-Higgs models [27, 28], Higgs portal models [29, 30], and models
inspired by warped extra dimensions [31]. In the last class of models, the new
heavy particles may have spin 0 (“radions”) or spin 2 (“gravitons”) [32]. In
this paper, the resulting “resonant” $\PH\PH$ production is sought for mass
values of X from 250 to 1000, and the width of X is assumed to be negligible
compared to the experimental resolution in $m_{\PH\PH}$. This creates a peak
in the reconstructed $m_{\PH\PH}$ distribution around the mass
$m_{\textrm{X}}$ of the resonance. The Feynman diagram for this process is
shown in Fig. 3. Resonance masses above 1are strongly constrained by searches
targeting boson decays to bottom quarks [33, 34], as the selection and
reconstruction efficiency, in particular of the trigger, for hadronic decays
increases and relevant backgrounds decrease with energy.

Figure 3:  Leading order Feynman diagram for resonant $\PH\PH$ production.

In this paper, we present a search for nonresonant as well as resonant
$\PH\PH$ production in final states with multiple reconstructed leptons, ,
electrons (), muons (), or hadronically decaying tau leptons (). The search is
based on LHC proton-proton ($\Pp\Pp$) collision data recorded by the CMS
experiment at a center-of-mass energy of 13, corresponding to an integrated
luminosity of 138. Signal candidate events are subdivided into seven mutually
exclusive “search categories” based on $\ell$ (, ) and multiplicity: two same-
sign $\ell$ with less than two ($2\ell\text{ss}$), three $\ell$ with no
($3\ell$), four $\ell$ ($4\ell$), three $\ell$ with one additional
($3\ell+1\tauh$), two $\ell$ with two ($2\ell+2\tauh$), one $\ell$ with three
($1\ell+3\tauh$), or four with no $\ell$ ($4\tauh$). In final states with a
total of four $\ell$ and , the charge sum of all $\ell$ and candidates is
required to be zero. The seven search categories target $\PH\PH$ signal events
in which the boson pair decays into $\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$, or
$\Pgt\Pgt\Pgt\Pgt$. Multivariate analysis (MVA) methods are used to
distinguish the $\PH\PH$ signal from backgrounds.

Phenomenological studies of the prospects for discovering $\PH\PH$ signal in
the $\PW\PW\PW\PW$ decay mode are documented in Refs. [35, 36, 37, 38, 39].
The ATLAS Collaboration published the results of a search for nonresonant as
well as resonant $\PH\PH$ production in this decay mode based on 36of $\Pp\Pp$
collision data recorded at $\sqrt{s}=13$ [40], placing an upper limit of 160
times the SM predicted cross section for nonresonant $\PH\PH$ production at
95% confidence level (). This paper presents the first search for boson pairs
decaying to $\PW\PW\Pgt\Pgt$ and $\Pgt\Pgt\Pgt\Pgt$. Searches for $\PH\PH$
production in $\Pp\Pp$ collisions at $\sqrt{s}=7$, 8, and 13have previously
been performed by the CMS and ATLAS Collaborations in the decay modes
$\PQb\PQb\PGg\PGg$ [41, 42], $\PQb\PQb\PQb\PQb$ [43, 34, 44],
$\PQb\PQb\PGt\PGt$ [45, 46], $\PQb\PQb\PW\PW$ [47, 48], and $\PW\PW\PGg\PGg$
[49]. Limits on $\PH\PH$ production obtained from a combination of some of
these analyses have been published by the CMS and ATLAS Collaborations [50,
51].

Searches targeting the $\PQb\PQb\PGg\PGg$ [41] and $\PQb\PQb\PQb\PQb$ [43]
final states in CMS, and $\PQb\PQb\PGg\PGg$ [42] in ATLAS, provide the
strongest constraints on nonresonant $\PH\PH$ production to date, with 95%
upper limits ranging from 3.9 to 7.7 times the SM predicted cross section. The
corresponding lower bounds on $\kappa_{\lambda}$ vary from $-3.3$ to $-1.5$,
with upper bounds between 6.7 and 9.4. The ATLAS $\PQb\PQb\PGg\PGg$ analysis
places a 95% upper limit of 0.64 pb on resonant $\PH\PH$ production with a
mass around 250, and lower limits for higher masses, down to 0.05 pb at 1. The
only published $\PH\PH$ search using an EFT approach comes from CMS in the
$\PQb\PQb\PGg\PGg$ final state, with 95% upper limits on the $\PH\PH$
production cross section ranging from 0.1 to 0.6 pb, depending on the EFT
scenario [41].

The paper is structured as follows. A brief overview of the CMS detector is
given in Section 0.2. In Section 0.3, we list the data sets and simulation
samples used. The reconstruction of , , , and jets, along with various
kinematic observables, is detailed in Section 0.4. This is followed by a
description, in Section 0.5, of the event selection criteria defining the
seven search categories. The multivariate methods used to distinguish the
$\PH\PH$ signal from backgrounds are detailed in Section 0.6. The estimation
of these backgrounds is described in Section 0.7, followed by an outline of
the relevant systematic uncertainties in Section 0.8. The statistical
procedure used to extract limits on the $\PH\PH$ production rate in the SM, as
well as constraints on SM coupling strengths, EFT benchmark scenarios, and
resonant $\PH\PH$ production rates are presented in Section 0.9. The paper
concludes with a summary in Section 0.10.

##  0.2 The CMS detector

The central feature of the CMS apparatus is a superconducting solenoid of 6m
internal diameter, providing a magnetic field of 3.8T. Within the solenoid
volume are a silicon pixel and strip tracker, a lead tungstate crystal
electromagnetic calorimeter (ECAL), and a brass and scintillator hadron
calorimeter (HCAL), each composed of one barrel and two endcap sections. The
silicon tracker measures charged particles within the pseudorapidity range
$\abs{\eta}<2.5$ for data recorded in 2016, and within the range
$\abs{\eta}<3.0$ for data recorded in 2017 and 2018. The ECAL is a fine-
grained hermetic calorimeter with quasi-projective geometry, and is divided
into a barrel region covering $\abs{\eta}<1.5$, and two endcaps that extend to
$\abs{\eta}=3.0$. The HCAL barrel and endcaps similarly cover the region
$\abs{\eta}<3.0$. Forward calorimeters extend beyond these endcaps to
$\abs{\eta}=5.0$. Muons are detected within the range $\abs{\eta}<2.4$ by gas-
ionization detectors embedded in the steel flux-return yoke outside the
solenoid. Collision events of interest are selected using a two-tiered trigger
system. The level-1 trigger, composed of custom hardware processors, uses
information from the calorimeters and muon detectors to select less than
100kHz of events from a 40MHz base event rate, within a fixed latency of 4
[52]. The second tier, known as the high-level trigger, is a processor farm
which runs a version of the full event reconstruction software optimized for
fast processing, and reduces the event rate to around 1kHz before data storage
[53]. A more detailed description of the CMS detector, together with a
definition of the coordinate system used and the most relevant kinematic
variables, can be found in Ref. [54].

##  0.3 Data samples and Monte Carlo simulation

The analyzed $\Pp\Pp$ collision data correspond to an integrated luminosity of
138, collected by the CMS detector over three years: 36in 2016, 42in 2017, and
60in 2018 [55, 56, 57]. Only the data-taking periods during which all detector
components were fully operational are included in the analysis. This analysis
uses triggers requiring one or more reconstructed , , or candidates to be
associated with the same collision vertex. The exact triggers and their
thresholds varied slightly from year to year because of changes in luminosity
and detector conditions, as well as improvements to the trigger algorithms.
The transverse momentum ($\pt$) thresholds imposed by the trigger on the
“leading” (highest ), “subleading” (second-highest ), and third , , or , and
the corresponding $\eta$ requirements for each year are shown in Table 0.3.
All triggers include identification and isolation requirements on the , , and
candidates [53]. When combined, the triggers achieve an efficiency of 95–100%
for simulated SM $\PH\PH$ signal events in each of the seven search
categories.

Selection requirements on and $\eta$ of reconstructed electrons (), muons (),
and hadronically decaying tau leptons () applied by the triggers used in this
analysis. The trigger thresholds for leading, subleading, and third , , or are
separated by commas. For trigger thresholds that varied over time, the range
of variation is indicated.  Trigger Selection requirements for reconstructed ,
, and objects Single  $\pt(\Pe)>27$–35 Single  $\pt(\PGm)>22$–27 Double
$\pt(\Pe)>23$, 12 \+  $\pt(\Pe)>23$, $\pt(\PGm)>8$ \+  $\pt(\PGm)>23$,
$\pt(\Pe)>8$–12 Double  $\pt(\PGm)>17$, 8 \+  $\pt(\Pe)>24$,
$\pt(\tauh)>20$–30, $\abs{\eta(\Pe,\tauh)}<2.1$ \+  $\pt(\PGm)>19$–20,
$\pt(\tauh)>20$–27, $\abs{\eta(\PGm,\tauh)}<2.1$ Double  $\pt(\tauh)>35$–40,
$\abs{\eta(\tauh)}<2.1$ Triple  $\pt(\Pe)>16$, 12, 8 Two +  $\pt(\Pe)>12$, 12,
$\pt(\PGm)>8$ Two +  $\pt(\PGm)>9$, 9, $\pt(\Pe)>9$ Triple  $\pt(\PGm)>12$,
10, 5

Monte Carlo (MC) simulated samples are used to model $\PH\PH$ signal events
and a wide range of SM background processes that produce final states with , ,
or . Background MC samples include processes producing a single or boson, two
bosons ($\PW\PW$, $\PW\PZ$, $\PZ\PZ$, $\PW\Pgg$, and $\PZ\Pgg$), three bosons
($\PW\PW\PW$, $\PW\PW\PZ$, $\PW\PZ\PZ$, $\PZ\PZ\PZ$, and $\PW\PZ\Pgg$), a
single boson (via gluon fusion, vector boson fusion, or associated production
with a $\PW$ or $\PZ$ boson), a single top quark, a top quark-antiquark pair
(), and top quarks associated with one or more bosons ($\PQt\PAQt\PW$,
$\PQt\PAQt\PZ$, $\PQt\PAQt\PH$, $\PQt\PH\PQq$, and $\PQt\PH\PW$). The $\PW\PW$
samples comprise same-sign and opposite-sign $\PW\PW$ pairs, including same-
sign $\PW\PW$ pairs produced by double-parton interactions. All samples that
include a boson were produced for a boson mass of 125. Most samples, including
the dominant $\PW\PZ$ and $\PZ\PZ$ backgrounds, were generated at next-to-LO
(NLO) and scaled to cross sections computed at NNLO in QCD.

A variety of $\PH\PH$ signal samples were generated at LO and at NLO accuracy
in QCD to simulate nonresonant $\PH\PH$ production, covering the
$\Pg\Pg\PH\PH$ and $\PQq\PQq\PH\PH$ production processes. Each of the bosons
was forced to decay to either $\PW\PW$, $\PZ\PZ$, or $\PGt\PGt$ in these
samples. The LO samples are used for the training of machine learning
algorithms, while the NLO samples are used to extract the rate of the $\PH\PH$
signal from the data. Separate $\Pg\Pg\PH\PH$ samples are produced for SM
$\PH\PH$ production and for a total of twelve benchmark (BM) scenarios in the
EFT approach [24]. These benchmarks, as well as the seven benchmarks from Ref.
[58], represent different combinations of $\kappa_{\lambda}$, $\kappa_{\PQt}$,
$\text{c}_{\Pg}$, $\text{c}_{2\Pg}$, and $\text{c}_{2}$ values, and are chosen
to probe distinct classes of $\PH\PH$ kinematic configurations. We will refer
to these benchmarks as JHEP04 BM1-12, and JHEP03 BM1-7, respectively. The
benchmark JHEP04 BM8 is complemented by a modified version of this benchmark,
published in Ref. [59], denoted as JHEP04 BM8a. The parameter values of these
twenty BM scenarios are shown in Table 0.3. The values of the $\text{c}_{\Pg}$
and $\text{c}_{2\Pg}$ couplings published in Ref. [58] have been scaled by
factors of 1.5 and $-3$, respectively, to convert them to the convention
introduced for these couplings in Ref. [24]. In order to increase the number
of simulated events and to model kinematic configurations not explicitly
generated, such as JHEP03 BM1-7, the $\Pg\Pg\PH\PH$ samples are merged and the
events in the merged samples are reweighted, using the procedure documented in
Ref. [60], to match the distributions in $m_{\PH\PH}$ and
$\abs{\cos\theta^{\ast}}$ computed at NLO accuracy and published in Ref. [59].
This procedure is applied to the LO and NLO $\Pg\Pg\PH\PH$ samples separately.
The symbol $\cos\theta^{\ast}$ denotes the cosine of the polar angle of one
with respect to the beam axis in the $\PH\PH$ rest frame. The $\PQq\PQq\PH\PH$
samples are produced only for SM $\PH\PH$ production.

Parameter values for $\kappa_{\lambda}$, $\kappa_{\PQt}$, $\text{c}_{2}$,
$\text{c}_{\Pg}$, and $\text{c}_{2\Pg}$ in MC samples modeling twenty
benchmark scenarios in the EFT approach, plus SM $\PH\PH$ production.
Benchmark $\kappa_{\lambda}$ $\kappa_{\PQt}$ $\text{c}_{2}$ $\text{c}_{\Pg}$
$\text{c}_{2\Pg}$ JHEP04 BM1 7.5 1.0 $-$1.0 0.0 0.0 JHEP04 BM2 1.0 1.0 0.5
$-$0.8 0.6 JHEP04 BM3 1.0 1.0 $-$1.5 0.0 $-$0.8 JHEP04 BM4 $-$3.5 1.5 $-$3.0
0.0 0.0 JHEP04 BM5 1.0 1.0 0.0 0.8 $-$1.0 JHEP04 BM6 2.4 1.0 0.0 0.2 $-$0.2
JHEP04 BM7 5.0 1.0 0.0 0.2 $-$0.2 JHEP04 BM8 15.0 1.0 0.0 $-$1.0 1.0 JHEP04
BM8a 1.0 1.0 0.5 4/15 0.0 JHEP04 BM9 1.0 1.0 1.0 $-$0.6 0.6 JHEP04 BM10 10.0
1.5 $-$1.0 0.0 0.0 JHEP04 BM11 2.4 1.0 0.0 1.0 $-$1.0 JHEP04 BM12 15.0 1.0 1.0
0.0 0.0 JHEP03 BM1 3.94 0.94 $-$1/3 0.75 $-$1 JHEP03 BM2 6.84 0.61 1/3 0 1
JHEP03 BM3 2.21 1.05 $-$1/3 0.75 $-$1.5 JHEP03 BM4 2.79 0.61 1/3 $-$0.75
$-$0.5 JHEP03 BM5 3.95 1.17 $-$1/3 0.25 1.5 JHEP03 BM6 5.68 0.83 1/3 $-$0.75
$-$1 JHEP03 BM7 $-$0.10 0.94 1 0.25 0.5 SM 1.0 1.0 0.0 0.0 0.0

Resonant $\PH\PH$ production was simulated at LO for both spin-0 (radion) and
spin-2 (graviton) scenarios with $m_{\textrm{X}}$ values of 250, 260, 270,
280, 300, 320, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, and
1000.

All MC samples were generated using either v2.2.2 or v2.3.3 [61, 62], v2 [63,
64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79], or v8.2 [80].
The production of $\PH\PH$ signal samples with is based on the work published
in Refs. [81, 82]. The parton distribution functions (PDFs) of the proton are
modeled using the NNPDF3.0 and NNPDF3.1 PDF sets [83, 84, 85, 86, 87]. Parton
shower, hadronization processes, and decays are modeled by , using the tunes
CP5, CUETP8M1, CUETP8M2, or CUETP8M2T4 [88, 89, 90], depending on the process
and the data-taking period that is being modeled. The matching of matrix
elements to parton showers is performed using the MLM scheme [91] for the LO
samples and the FxFx scheme [92] for the NLO samples. The interactions of
particles with the CMS detector material was simulated in detail using [93].
Simulated events were reconstructed using the same procedure as in data. The
response of the trigger is included in the simulation. Additional $\Pp\Pp$
interactions (pileup) were generated with and overlaid on all MC events, with
event weights used to match the collision multiplicity to the distribution
inferred from data. Residual differences between data and simulation are
rectified by applying corrections to simulated events.

##  0.4 Event reconstruction

The CMS particle-flow (PF) algorithm [94] aims to reconstruct and identify
each individual particle in an event, using an optimized combination of
information from the various elements of the CMS detector. The particles are
subsequently classified into five mutually exclusive types: electrons, muons,
photons, and charged and neutral hadrons. These particles are then combined to
reconstruct hadronic decays, jets, and the missing transverse momentum in the
event.

The candidate vertex with the largest value of summed physics-object $\pt^{2}$
is taken to be the primary $\Pp\Pp$ interaction vertex. The physics objects
used for this determination are the jets, clustered using the infrared and
collinear safe anti-algorithm [95, 96], with the tracks assigned to candidate
vertices as inputs, and the associated missing transverse momentum, taken as
the negative vector sum of the of those jets.

Electrons are reconstructed within the geometric acceptance of the tracking
detectors ($\abs{\eta}<2.5$) by combining information from the tracker and the
ECAL [97]. They are initially identified using an MVA classifier which
distinguishes real electrons from hadrons, along with requirements that the
track be associated with the collision vertex, and limits on hadronic energy
deposits separated by $\Delta R<0.4$ from the electrons (their “isolation”).
The angular separation between two particles is defined as $\Delta
R=\sqrt{(\eta_{1}-\eta_{2})^{2}+(\phi_{1}-\phi_{2})^{2}}$, where the symbol
$\phi$ refers to the azimuthal direction of the particle. Electrons passing
this initial selection are referred to as “loose”. In this analysis, events
with electrons originating from hadron decays (“nonprompt”), or with hadrons
misidentified as electrons, constitute the largest source of background. This
motivates the use of an additional MVA classifier, which is trained to select
“prompt” electrons from , , and lepton decays, and to reject nonprompt or
misidentified electrons. This MVA classifier was previously used for
measurements of $\PQt\PAQt\PH$ production in events with multiple leptons
[98]. It combines observables comparing measurements of the electron energy
and direction in the tracker and the ECAL, the compactness of the electron
cluster, the bremsstrahlung emitted along the electron trajectory, and the
electron isolation. Two levels of thresholds on the output of this MVA
classifier are used in the analysis, referred to as the “tight” and “medium”
electron selections for the more and less restrictive thresholds,
respectively. The tight selection has an average efficiency of 60% for
electrons from SM $\PH\PH$ decays. Only the electrons passing the tight
selections are used to reconstruct signal candidate events, while data events
passing the medium electron selections and failing the tight selections are
used to estimate the contribution of misidentified- and nonprompt-electron
backgrounds in each search category. Compared to Ref. [98], this analysis uses
lower thresholds on the MVA classifier output for the medium and tight
electron selections, in order to increase the efficiency in particular for
low-electrons, which frequently appear in the $\PH\PH$ signal events studied
in this analysis. Electrons from photon conversions in the tracker are
suppressed by requiring that the track is missing no hits in the innermost
layers of the silicon tracker, and is not matched to a reconstructed
conversion vertex. In the $2\ell\text{ss}$ category, further electron
selection criteria are applied, which require agreement among three
independent measurements of the electron charge, including the Gaussian sum
filter and Kalman filter track curvatures, as well as the ECAL supercluster
position [99]. The remaining charge misidentification rate is measured to be
less than 0.1% for $\abs{\eta}<1.479$, and under 0.4% for $\abs{\eta}>1.479$.
The charge quality requirement reduces the electron identification efficiency
by about $4\%$.

Muons are reconstructed by extrapolating tracks in the silicon tracker to hits
in the gas-ionization muon detectors embedded in the steel flux-return yoke
outside the solenoid [100]. To pass the initial loose identification
requirement for this analysis, muons must satisfy criteria related to
isolation and track proximity to the primary interaction vertex, as well as
track quality observables and matching between the tracker and muon chambers.
Additional requirements on the prompt vs. nonprompt muon identification MVA
classifier from Ref. [98] serve to select muons passing a tight selection for
signal candidate events, and a medium selection for nonprompt background
estimation. Inputs to this MVA classifier include energy deposits close to the
muon in the ECAL and HCAL, the hits and track segments reconstructed in the
muon detectors located outside the solenoid, the quality of the spatial
matching between the track segments reconstructed in the silicon tracker and
in the muon detectors, and the isolation of the muon with respect to other
particles. Again, lower selection thresholds on the MVA classifier output
compared to Ref. [98] bring higher efficiency for the $\PH\PH$ signal,
amounting to 80% per muon in simulated SM $\PH\PH$ events for the tight
selection. In the $2\ell\text{ss}$ channel, the uncertainty in the curvature
of the muon track is required to be less than 20% to ensure a high-quality
charge measurement [98]. This requirement reduces the muon identification
efficiency by about 2%.

Hadronic decays of tau leptons are identified using the “hadrons-plus-strips”
algorithm [101]. This algorithm classifies individual hadronic decay modes of
the by combining charged hadrons from the PF reconstruction with neutral
pions. The latter are reconstructed by clustering electrons and photons into
rectangular strips, which are narrow in $\eta$ but wide in the $\phi$
direction. The spread in $\phi$ accounts for photons originating from neutral
pion decays that convert into electron-positron ($\Pem\Pep$) pairs while
traversing the silicon tracker. The $\Pem$ and $\Pep$ are bent in opposite
directions in $\phi$ by the magnetic field, and may further emit
bremsstrahlung photons before reaching the ECAL. The decay modes considered in
this analysis produce one charged pion or kaon plus zero to two neutral pions
(collectively referred to as “one-prong” ), or three charged pions or kaons
plus zero or one neutral pion (referred to as “three-prong” ). The DeepTau
algorithm [102] distinguishes true objects from quark and gluon jets,
electrons, and muons using a convolutional artificial neural network (ANN)
[103] with 42 high-level observables as input, together with low-level
information obtained from the silicon tracker, ECAL, HCAL, and the muon
detectors. The former include the , $\eta$, $\phi$, and mass of the candidate,
the reconstructed decay mode, its isolation with respect to charged and
neutral particles, and the estimated distance that the lepton traverses
between its production and decay. For three-prong candidates, this distance is
determined by reconstructing the decay vertex, while for one-prong candidates,
the transverse impact parameter of the charged pion track with respect to the
primary $\Pp\Pp$ interaction vertex is used as an estimate of the distance.
The low-level information quantifies the particle activity within two
$\eta\times\phi$ grids, centered on the direction of the candidate: an inner
grid of size $0.2\times 0.2$, filled with $0.02\times 0.02$ cells, and an
outer grid of size $0.5\times 0.5$ (partially overlapping with the inner
grid), with $0.05\times 0.05$ cells. Selected candidates in this analysis must
have $\pt>20$and $\abs{\eta}<2.3$, and are subjected to two levels of
thresholds on the output of the ANN that separates from quark and gluon jets,
referred to as the tight and medium selections, respectively.

Hadronic jets (j) are reconstructed with the anti-algorithm using the
particles reconstructed with the PF algorithm as input, and serve to identify
$\PH\to\PW\PW\to\text{j}\text{j}\ell\PGn$ decays in this analysis. Jets
reconstructed with distance parameters of 0.4 (“small-radius jets”) and 0.8
(“large-radius jets”) are both used: two small-radius jets to reconstruct the
two quarks from low-boson decays, or a single large-radius jet to reconstruct
high-boson decays, where the quarks are collimated.

The effect of pileup on the reconstruction of large-radius jets is mitigated
by applying the pileup per particle identification algorithm (PUPPI) [104,
105] to the collection of particles used as input to the jet reconstruction.
Electrons and muons passing the loose selection are also removed from the
input collection, so that leptons produced in
$\PH\to\PW\PW\to\to\text{j}\text{j}\ell\PGn$ decays of Lorentz-boosted bosons
are not clustered into large-radius jets. For small-radius jets, the effect of
pileup is reduced by removing charged particles identified with pileup
vertices from the jet reconstruction, and applying corrections to the jet
energy to account for neutral particles from pileup. Small-radius jets which
contain at least one PF particle matched to an electron, muon, or constituent
of a passing the medium selection criteria are discarded.

After calibration, the jet energy resolution at the central rapidities amounts
to 15–20% at 30, 10% at 100, and 5% at 1 [106]. This analysis only considers
jets reconstructed in the region $\abs{\eta}<2.4$. Small-radius jets must have
$\pt>25$, while large-radius jets must have $\pt>170$. Additional criteria
requiring that each large-radius jet contain exactly two identifiable,
energetic subjets are applied to specifically select those from boosted
hadronic boson decays [107].

Events containing small-radius jets identified with the hadronization of
bottom quarks ($\PQb$ jets) are vetoed in this analysis. The DeepJet algorithm
[108] exploits observables related to the long lifetime of $\PQb$ hadrons and
the higher particle multiplicity and mass of $\PQb$ jets compared to light
quark and gluon jets. Both “loose” and “medium” $\PQb$ jet selections on the
DeepJet output are employed in this analysis, corresponding to $\PQb$ jet
selection efficiencies of 84 and 70%, while the misidentification rates for
light-quark or gluon jets are 11 and 1.1%, respectively.

The missing transverse momentum vector $\ptvecmiss$ is computed as the
negative vector $\pt$ sum of all the particles reconstructed by the PF
algorithm in an event, and its magnitude is denoted as $\ptmiss$ [109]. The
$\ptvecmiss$ is modified to account for corrections to the energy scale of the
reconstructed jets in the event. A linear discriminant, denoted as
$\pt^{\text{miss,LD}}$, is employed to remove background events in which the
reconstructed arises from resolution effects. The discriminant is defined by
the relation $\pt^{\text{miss,LD}}=0.6\ptmiss+0.4\HT^{\text{miss}}$, where
$\HT^{\text{miss}}$ corresponds to the magnitude of the vector sum of , , and
passing the medium selection criteria, and small-radius jets satisfying the
criteria detailed above [110].

##  0.5 Event selection

Events are selected with the aim of maximizing the acceptance for $\PH\PH$
decays to $\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$, and $\Pgt\Pgt\Pgt\Pgt$, while
simultaneously rejecting the large backgrounds from multijet production,
single and pair production of and bosons, and production. To achieve this, we
require multiple reconstructed $\ell$ or to be associated with the primary
interaction vertex. The $\ell$ and may originate from the decay of a $\PW$
boson or a $\PGt$ lepton. Seven mutually exclusive search categories,
distinguished by the number of reconstructed $\ell$ and candidates, are
included in the analysis: $2\ell\text{ss}$, $3\ell$, $4\ell$, $3\ell+1\tauh$,
$2\ell+2\tauh$, $1\ell+3\tauh$, and $4\tauh$. Here “ss” indicates a same-sign
$\ell\ell$ pair, with two leptons of identical electric charge. The $\ell$ and
candidates selected in any of the seven search categories must pass the tight
selection criteria described in Section 0.4. In addition, they are required to
pass category-specific thresholds motivated by the trigger selection. Further
requirements are placed on the sum of $\ell$ and charges, and, in two
categories, on the discriminant $\pt^{\text{miss,LD}}$ and the multiplicity of
jets.

The leading and subleading leptons in the $2\ell\text{ss}$ category must pass
selection thresholds of 25 and 15, respectively. Events in this category are
required to contain two or more small-radius jets, or at least one large-
radius jet, targeting hadronic boson decays. Dielectron events must have
$\pt^{\text{miss,LD}}>30$and $m(\ell\ell)<81$or $m(\ell\ell)>101$, in order to
suppress charge-misidentified $\PZ\to\Pe\Pe$ background. If the event contains
a , the charge of the must be opposite to the charge of the leptons. After
this selection, the main backgrounds in the $2\ell\text{ss}$ category arise
from $\PW\PZ$ production, from $\PW\Pgg$ events in which the photon converts
into an $\Pem\Pep$ pair and either the $\Pem$ or the $\Pep$ fails to get
reconstructed, and from events in which one or both reconstructed leptons are
due to a nonprompt $\ell$ or misidentified hadron, as shown in Table 0.5. The
“other” background given in the table is dominated by same-sign boson pairs
and $\PW\PW\PW$ production. The $\PW\PW\PW\PW$ decay mode accounts for roughly
70% of SM $\PH\PH$ signal events selected in the $2\ell\text{ss}$ category,
with $\PW\PW\Pgt\Pgt$ events accounting for the other 30%.

In the $3\ell$ category, the leading, subleading, and third $\ell$ are
required to have values greater than 25, 15, and 10, respectively, and the sum
of their charges must be either $+1$ or $-1$. At least one small- or large-
radius jet must be present, and the $\pt^{\text{miss,LD}}$ quantity must be
greater than 30, or 45if there is at least one SFOS $\ell\ell$ pair in the
event. Again, backgrounds are dominated by $\PW\PZ$ production and events with
misidentified $\ell$. Notable contributions to the “other” background arise
from $\PW\PW\PW$ and $\PW\PW\PZ$ production. The signal composition is similar
to the $2\ell\text{ss}$ category.

The $4\ell$ category has identical lepton selection criteria to the $3\ell$
category, except that the third $\ell$ must have $\pt>15$, and a fourth $\ell$
with $\pt>10$is required, and the sum of the four lepton charges is required
to be equal to zero. In this category and all the remaining categories, there
are no selection requirements on jets or $\pt^{\text{miss,LD}}$. Almost $70\%$
of signal events come from the $\PW\PW\PW\PW$ decay mode, and about 30% from
$\PW\PW\Pgt\Pgt$, while $\PZ\PZ$ production accounts for 85% of the
background.

Events in the $3\ell+1\tauh$ category are required to satisfy the $3\ell$
criteria on the $\ell$ objects, except that an additional with $\pt>20$and
charge opposite to the sum of the $\ell$ charges is required. Background
events in which the reconstructed fails a loose selection on the output of the
ANN of the DeepTau algorithm that separates from electrons, or falls near the
ECAL barrel-endcap transition region in $1.460<\abs{\eta}<1.558$ are removed.
About 70% of signal events come from the $\PW\PW\Pgt\Pgt$ decay mode, while
$\PZ\PZ$ production and events with at least one misidentified $\ell$ or
dominate the background.

In the $2\ell+2\tauh$ category, the leading and subleading $\ell$ are required
to pass thresholds of 25 and 15, while the two must have $\pt>20$. The sum of
$\ell$ plus charges is required to be zero. Signal contributions are mostly
from the $\PW\PW\Pgt\Pgt$ (60%) and $\Pgt\Pgt\Pgt\Pgt$ (40%) decay modes,
while background contributions arise from $\PZ\PZ$ production and events with
a misidentified $\ell$ or candidate.

In the $1\ell+3\tauh$ category, the leading $\ell$ is required to satisfy the
conditions $\abs{\eta}<2.1$ and $\pt>20$ (15)if it is an electron (muon). The
leading, subleading, and third must have $\pt>40$, 30, and 20, respectively,
and the sum of and $\ell$ charges is required to be zero. Background events
containing a $\PZ\to\Pe\Pe$ decay where one electron is misidentified as a are
vetoed by discarding events containing an -pair of opposite charge and mass
$71<m(\Pe\tauh)<101$, and in which the either fails a loose selection on the
discriminant that separates from electrons, or falls into the region
$1.460<\abs{\eta}<1.558$. Around 80% of $\PH\PH$ signal events selected in the
$1\ell+3\tauh$ category are from $\Pgt\Pgt\Pgt\Pgt$ and 20% from the
$\PW\PW\Pgt\Pgt$ decay mode, while the majority of background events stem from
$\PZ\PZ$ production or contain a misidentified $\ell$ or .

The $4\tauh$ category has identical requirements on the leading three , and
the fourth is also required to have $\pt>20$. Given the extremely low
backgrounds in this category, no charge sum criterion or $\PZ\to\Pe\Pe$ veto
is applied. Almost all signal events come from the $\Pgt\Pgt\Pgt\Pgt$ decay
mode, while 55% of the background events contain at least one misidentified
candidate, and the remainder arises from $\PZ\PZ$ (30%) and single Higgs boson
($15\%$) production.

In all seven search categories, the background contamination from processes
with top quarks is reduced by discarding events with at least one selected
small-radius jet passing the medium $\PQb$ jet identification, or at least two
passing the loose $\PQb$ jet identification. Leptons originating from low-mass
Drell–Yan production, decays of $\JPsi$ and $\Upsilon$ mesons, cascade decays
of bottom quarks, and photon conversions are removed by vetoing events
containing any pair of loose $\ell$ with mass $m(\ell\ell)<12$. To eliminate
overlap with events selected in the ongoing search for $\PH\PH$ production in
the $\PQb\PAQb\PZ\PZ$, $\PZ\PZ\to 4\ell$ decay mode, no event in the
$2\ell\text{ss}$, $3\ell$, and $4\ell$ categories may contain two same-flavor
opposite-sign (SFOS) loose $\ell\ell$ pairs with a combined mass less than
140. In addition, to reduce the $\PZ\to\ell\ell$ background, these three
categories along with $2\ell+2\tauh$ and $3\ell+1\tauh$ exclude events where
any SFOS loose $\ell\ell$ pair has an invariant mass of 81–101(boson veto).

A summary of the event selection criteria applied in the different categories
is given in Table 0.5. Criteria that are common to all seven search categories
are given in Table 0.5.

The number of events selected in the signal regions of each of the seven
search categories, and in two control regions (CRs) which are used to validate
the modeling of the $\PW\PZ$ and $\PZ\PZ$ backgrounds, are given in Table 0.5.
The contribution expected from nonresonant $\PH\PH$ production with event
kinematics as predicted by the SM, but 30 times the SM cross section, is given
separately for $\PH\PH$ decays into $\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$, and
$\Pgt\Pgt\Pgt\Pgt$ in the upper three rows of each table. The event yields
given in the rows labeled $\PW\PW\PW\PW$ include a small contribution from
$\PH\PH$ decays into $\PW\PW\PZ\PZ$ and $\PZ\PZ\PZ\PZ$, and, similarly, the
numbers quoted in the rows labeled $\PW\PW\Pgt\Pgt$ include a small
contribution from $\PH\PH$ decays into $\PZ\PZ\Pgt\Pgt$. The $3\ell~{}\PW\PZ$
and $4\ell~{}\PZ\PZ$ CRs are described in Section 0.7.

Reconstructed object and event selection requirements in all seven search
categories. Electrons or muons in the $\ell\ell$ pairs include any leptons
passing the loose selection criteria.  Object and event properties Selection
criteria Lepton and pseudorapidity $\abs{\eta}<2.5$ for , $\abs{\eta}<2.4$ for
, $\abs{\eta}<2.3$ for  Dilepton invariant mass $m_{\ell\ell}>12$(all
$\ell\ell$ pairs) Four-lepton invariant mass $m_{4\ell}>140$(any two SFOS
$\ell\ell$ pairs) $\PQb$ jet veto 0 medium and $\leq 1$ loose $\PQb$-tagged
small-radius jet

Event selection criteria applied in the seven search categories. The
thresholds for $\ell$ and with the highest, second-, third-, and fourth-
highest are separated by slashes. The symbol “” indicates that no requirement
is applied.  Category $2\ell\text{ss}$ $3\ell$ $4\ell$ Targeted $\PH\PH$
decays $\PW\PW\PW\PW$ $\PW\PW\PW\PW$ $\PW\PW\PW\PW$ Trigger Single- and
Single-, double- Single-, double- double-lepton and triple-lepton and triple-
lepton Lepton  $>$25 / 15 $>$25 / 15 / 10 $>$25 / 15 / 15 / 10 Lepton charge
sum $\pm 2$, with charge quality $\pm 1$ 0 requirements applied Dilepton
invariant mass $\abs{m_{\ell\ell}-m_{\PZ}}>10$†
$\abs{m_{\ell\ell}-m_{\PZ}}>10$‡ $\abs{m_{\ell\ell}-m_{\PZ}}>10$‡ Jets $\geq$2
small-radius jets or $\geq$1 small-radius jet or $\geq$1 large-radius jet
$\geq$1 large-radius jet Missing $\pt$ $\pt^{\text{miss,LD}}>30$§
$\pt^{\text{miss,LD}}>30$∥ Category $3\ell+1\tauh$ $2\ell+2\tauh$ Targeted
$\PH\PH$ decays $\PW\PW\Pgt\Pgt$ $\PW\PW\Pgt\Pgt$, $\Pgt\Pgt\Pgt\Pgt$ Trigger
Single-, double-, Single- and and triple-lepton double-lepton Lepton  $>$25 /
15 / 10 $>$25 / 15 $>$20 $>$20 Lepton and charge $\ell$ and charges sum to 0
$\ell$ and charges sum to 0 Dilepton invariant mass
$\abs{m_{\ell\ell}-m_{\PZ}}>10$‡ $\abs{m_{\ell\ell}-m_{\PZ}}>10$‡ Category
$1\ell+3\tauh$ $4\tauh$ Targeted $\PH\PH$ decays $\Pgt\Pgt\Pgt\Pgt$
$\Pgt\Pgt\Pgt\Pgt$ Trigger Single-lepton, lepton+ Double- and double- Lepton
$\eta$ $\abs{\eta}<2.1$ Lepton  $>$20() or $>$15() $>$40 / 30 / 20 $>$40 / 30
/ 20 / 20 Lepton and charge $\ell$ and charges sum to 0 charges sum to 0
$\PZ\to\Pe\Pe$ veto $\abs{m_{\Pe\tauh}-86\GeV}>15$¶ † Applied to all SFOS
$\ell\ell$ pairs and electron pairs with the same charge.
‡ Applied to all SFOS $\ell\ell$ pairs.
§ Only applied to events containing two electrons.
∥ Tightened to $\pt^{\text{miss,LD}}>45$if event contains a SFOS $\ell\ell$
pair.
¶ For classified as electrons by the DeepTau algorithm or with
$1.460<\abs{\eta}<1.558$.

The number of expected and observed events in each of the seven search
categories, and in two control regions (CRs), which validate the modeling of
the $\PW\PZ$ and $\PZ\PZ$ backgrounds. The $\PH\PH$ signal represents the sum
of the $\Pg\Pg\PH\PH$ and $\PQq\PQq\PH\PH$ production processes and is
normalized to 30 times the event yield expected in the SM, corresponding to a
cross section of about 1 pb. The expected event yields are computed for the
values of nuisance parameters obtained from the ML fit described in Section
0.9. Quoted uncertainties represent the sum of statistical and systematic
components. Uncertainties that are smaller than half the value of the least
significant digit have been rounded to zero.  Process $2\ell\text{ss}$ $3\ell$
$4\ell$ SM $\PH\PH\to\PW\PW\PW\PW$ ($\times\,30$) 73 $\pm$ 2 33 $\pm$ 1 2.2
$\pm$ 0.1 SM $\PH\PH\to\PW\PW\Pgt\Pgt$ ($\times\,30$) 31 $\pm$ 1 12 $\pm$ 0
1.0 $\pm$ 0.1 SM $\PH\PH\to\Pgt\Pgt\Pgt\Pgt$ ($\times\,30$) 3 $\pm$ 0 1 $\pm$
0 0.1 $\pm$ 0.0 $\PW\PZ$ 2003 $\pm$ 19 1321 $\pm$ 13 0.4 $\pm$ 0.1 $\PZ\PZ$
121 $\pm$ 1 109 $\pm$ 1 54.7 $\pm$ 1.6 Misidentified $\ell$ 3939 $\pm$ 83 670
$\pm$ 26 2.3 $\pm$ 1.0 Conversion electrons 1009 $\pm$ 66 146 $\pm$ 12 0.9
$\pm$ 0.4 Electron charge misid. 366 $\pm$ 17 $<$1 $<$0.1 Single Higgs boson
216 $\pm$ c2 62 $\pm$ 1 2.4 $\pm$ 0.3 Other backgrounds 2592 $\pm$ 79 249
$\pm$ 11 1.1 $\pm$ 0.2 Total expected background 10 346 $\pm$ 91 2601 $\pm$ 32
64.8 $\pm$ 1.8 Data 10 344 2621 62 Process $3\ell+1\tauh$ $2\ell+2\tauh$
$1\ell+3\tauh$ $4\tauh$ SM $\PH\PH\to\PW\PW\PW\PW$ ($\times\,30$) 0.9 $\pm$
0.1 0.2 $\pm$ 0.0 0.2 $\pm$ 0.0 0.3 $\pm$ 0.0 SM $\PH\PH\to\PW\PW\Pgt\Pgt$
($\times\,30$) 4.1 $\pm$ 0.2 3.9 $\pm$ 0.4 0.6 $\pm$ 0.1 0.1 $\pm$ 0.0 SM
$\PH\PH\to\Pgt\Pgt\Pgt\Pgt$ ($\times\,30$) 0.9 $\pm$ 0.0 2.3 $\pm$ 0.1 2.6
$\pm$ 0.3 1.3 $\pm$ 0.2 $\PW\PZ$ 0.2 $\pm$ 0.0 $<$0.1 $<$0.1 $<$0.1 $\PZ\PZ$
24.3 $\pm$ 0.3 18.5 $\pm$ 0.4 1.9 $\pm$ 0.1 0.7 $\pm$ 0.1 Misidentified $\ell$
and  25.1 $\pm$ 3.1 33.5 $\pm$ 3.8 2.1 $\pm$ 1.6 1.5 $\pm$ 0.9 Conversion
electrons 0.1 $\pm$ 0.0 0.1 $\pm$ 0.1 $<$0.1 $<$0.1 Single Higgs boson 3.8
$\pm$ 0.2 2.9 $\pm$ 0.5 0.8 $\pm$ 0.3 0.4 $\pm$ 0.1 Other backgrounds 1.1
$\pm$ 0.1 1.2 $\pm$ 0.2 $<$0.1 $<$0.1 Total expected background 56.2 $\pm$ 3.3
57.0 $\pm$ 3.8 4.9 $\pm$ 1.6 2.6 $\pm$ 0.9 Data 55 55 $6$ $1$ Process
$3\ell~{}\PW\PZ$ CR $4\ell~{}\PZ\PZ$ CR $\PW\PZ$ 12 546 $\pm$ 55 $<$1 $\PZ\PZ$
799 $\pm$ 7 2032 $\pm$ 20 Misidentified $\ell$ 908 $\pm$ 41 13 $\pm$ 4
Conversion electrons 134 $\pm$ 9 3 $\pm$ 0 Other backgrounds 505 $\pm$ 17 44
$\pm$ 2 Total expected background 15 006 $\pm$ 64 2108 $\pm$ 14 Data 14 994
2096

##  0.6 Analysis strategy

The rate of the $\PH\PH$ signal is extracted through a binned maximum
likelihood (ML) fit to the distributions in the output of boosted decision
tree (BDT) classifiers [111], which are trained to discriminate the $\PH\PH$
signal from backgrounds. The data from each of the three years are fit
separately. Three classifiers are trained for each of the seven search
categories using a mix of MC simulation from all three years, targeting
nonresonant $\PH\PH$ production and resonant $\PH\PH$ production from the
decay of heavy particles of spin 0 and of spin 2. In the two categories with
high event yields ($2\ell\text{ss}$ and $3\ell$), and in the $4\ell$ category,
the BDT output binning is chosen such that each bin contains the same number
of expected $\PH\PH$ signal events. The number of bins is determined by the
condition that the relative statistical uncertainty in the background
prediction in each bin does not exceed 15%. In the remaining four categories
($3\ell+1\tauh$, $2\ell+2\tauh$, $1\ell+3\tauh$, and $4\tauh$), which have low
event yields, the binning for each BDT output distribution is chosen such that
a similar number of expected background events is contained in each bin.
Higher bin numbers correspond to a higher BDT output value, and feature a
higher signal-to-background ratio. For the SM $\PH\PH$ signal, the bins with
the highest BDT output values feature a signal-to-background ratio up to 10
times higher than the inclusive ratio in each category.

The inputs to the BDT classifiers include the and $\eta$ of reconstructed
$\ell$ and ; the angular separation $\Delta R$ within $\ell\ell$, $\ell\tauh$,
and $\tauh\tauh$ pairs; the $\Delta R$ between an $\ell$ or candidate and the
nearest jet; the number of jets in the event; the discriminant
$\pt^{\text{miss,LD}}$; the scalar sum of all reconstructed , , , and jets;
the “visible” mass of the Higgs boson pair, given by the mass of the system of
reconstructed , , , and jets; and the “full” mass of the $\PH\PH$ system,
including neutrinos, reconstructed using the algorithm described in Ref.
[112]. This algorithm targets $\PH\PH$ signal events decaying to
$\Pgt\Pgt\Pgt\Pgt$ and thus works best in the $4\tauh$ and $1\ell+3\tauh$
search categories. Distributions in some of the observables used as inputs to
the BDT classifiers in the $2\ell\text{ss}$ and $3\ell$ categories are shown
in Fig. 4.

Figure 4:  Distributions in a few observables used as inputs to the BDT
classifiers in the $2\ell\text{ss}$ and $3\ell$ categories: the scalar sum,
denoted as , of the two reconstructed $\ell$ and all small-radius jets in the
$2\ell\text{ss}$ category (upper left); the angular separation $\Delta R$
between the two $\ell$ in the $2\ell\text{ss}$ category (upper right); the
angular separation between $\ell_{3}$ and the nearest small-radius jet in the
$3\ell$ category (lower left); and $\pt^{\text{miss,LD}}$ in the $3\ell$
category (lower right). The $\ell_{3}$ in the $3\ell$ category is defined as
the $\ell$ that is not part of the opposite-sign $\ell\ell$ pair of lowest
mass. The distributions expected for the different background processes are
shown for the values of nuisance parameters obtained from the background-only
ML fit, in which the $\PH\PH$ signal is constrained to be zero.

These observables are complemented by further inputs, which allow the BDT to
learn that the distributions in the observables for $\PH\PH$ signal events
change as a function of the model parameters: the Higgs boson couplings
$\lambda$, $\text{y}_{\PQt}$, $\text{c}_{\Pg}$, $\text{c}_{2\Pg}$, and
$\text{c}_{2}$ for nonresonant $\PH\PH$ production, and the mass of the heavy
particle X in resonant $\PH\PH$ production. The training is performed using
simulated samples of signal and background events. Background contributions
arising from the misidentification of $\ell$ and candidates and from the
mismeasurement of the electron charge are included in the simulation. When
training the BDT that targets nonresonant $\PH\PH$ production, the values for
the couplings are chosen according to the twelve EFT benchmark scenarios given
in Ref. [24]. The BDT for nonresonant $\PH\PH$ production is trained on
$\PH\PH$ samples corresponding to the SM prediction as well as these twelve
benchmark scenarios, indicated by thirteen binary inputs to the BDT. The BDT
classifiers used for the analysis of resonant $\PH\PH$ production are trained
on the full set of resonance masses listed in Section 0.3, and the resonance
mass is used as an input to the BDT.

Each simulated background event is replicated thirteen times in the training
sample when training the BDT for nonresonant $\PH\PH$ production, and nineteen
times when training the BDTs that target resonant $\PH\PH$ production, setting
the binary inputs for the Higgs boson couplings and the mass of the heavy
particle X to different values in each duplicate. The number of training
events is increased by applying the medium $\ell$ and identification criteria
instead of the tight ones. Weights are applied to background events arising
from different sources, such that the relative fractions of different types of
backgrounds in the training match the fractions expected in the signal region
of the analysis, when the tight $\ell$ and identification criteria are
applied. The MC samples used for the BDT training overlap with the samples
used to model signal and background contributions in the analysis. To avoid
potential biases, we split the training samples into two samples of equal
size, based on even and odd event numbers. The BDTs trained on even events are
evaluated on odd events, and vice versa, thereby ensuring that BDT are not
trained and evaluated on the same events. The training is performed using the
XGBoost algorithm [113], interfaced to the Scikit-learn machine learning
library [114]. The parameters of the BDT training (so-called
“hyperparameters”) are optimized using the particle swarm optimization
algorithm described in Ref. [115].

##  0.7 Background estimation

Background contributions are classified as either “reducible” or
“irreducible”. In this analysis, three types of reducible backgrounds are
considered, arising from misidentified $\ell$ or , electron charge
misidentification (“flips”), and electrons from photon conversions. Background
events in which all selected $\ell$ and come from , , or boson decays, and are
reconstructed with the correct charge, are considered “irreducible”. The
$\ell$/misidentification and electron charge misidentification backgrounds are
both determined from data, while electron conversions and irreducible
backgrounds are modeled using MC simulation.

The $\ell$/misidentification background (which includes nonprompt leptons) is
the largest reducible background in all search categories. Nonprompt $\ell$
are either electrons or muons produced in bottom and charm quark decays, or
muons that originate from pion and kaon decays. Hadron jets may also be
misidentified as electrons or . The $\ell$/misidentification background
estimate is detailed in Section 0.7.1. The electron charge misidentification
background is only relevant for the $2\ell\text{ss}$ search category, and is
described in Section 0.7.2. The modeling of photon conversion events by the MC
simulation has been validated in data as described in Ref. [98].

The main contribution to the irreducible background arises from $\PW\PZ$
production in the $2\ell\text{ss}$ and $3\ell$ categories, and $\PZ\PZ$
production in the remaining five categories. The production of pairs of bosons
(, , , or ) other than $\PW\PZ$, $\PZ\PZ$ and $\PH\PH$, and production of
bosons with top quarks, including $\PW\PGg$, $\PZ\PGg$, $\PW\PH$, $\PZ\PH$,
$\PQt\PH$, $\PQt\PAQt\PH$, $\PQt\PW$, $\PQt\PAQt\PW$, $\PQt\PZ$,
$\PQt\PAQt\PZ$, $\PQt\PGg$, and $\PQt\PAQt\Pgg$, constitute subdominant
additional backgrounds. The $\PQt\PZ$ and $\PQt\PAQt\PZ$ backgrounds also
include contributions from off-shell $\PQt\PAQt\PGg^{\ast}$ and
$\PQt\PGg^{\ast}$ production. Background processes which include at least one
top quark are suppressed by the $\PQb$ jet veto described in Section 0.5, but
are still sizable compared to the expected $\PH\PH$ signal. All irreducible
backgrounds are modeled using the MC simulation.

The modeling of the dominant irreducible $\PW\PZ$ and $\PZ\PZ$ backgrounds is
validated in dedicated CRs, the “$3\ell~{}\PW\PZ$” CR and “$4\ell~{}\PZ\PZ$”
CR. These CRs are based on the signal regions of the $3\ell$ and $4\ell$
categories, but with the boson veto inverted. The $3\ell~{}\PW\PZ$ and
$4\ell~{}\PZ\PZ$ CRs are included in the ML fit that is used to extract the
$\PH\PH$ signal, described in Section 0.9, thereby providing in-situ
constraints on the $\PW\PZ$ and $\PZ\PZ$ backgrounds and on systematic
uncertainties related to lepton identification and trigger efficiency.
Distributions in kinematic observables in the $3\ell~{}\PW\PZ$ and
$4\ell~{}\PZ\PZ$ CRs are shown in Fig. 5. The transverse mass,
$\mT=\sqrt{2\,\pt^{\ell}\,\ptmiss\left(1-\cos\Delta\phi\right)}$, in the
$3\ell~{}\PW\PZ$ CR is computed using the $\ell$ that is not identified as
originating from the boson decay. The symbol $\Delta\phi$ refers to the angle
in the transverse plane between the $\ell$ momentum and the . The observable
$m_{4\ell}$ refers to the mass of the $4\ell$ system in the $4\ell~{}\PZ\PZ$
CR.

Figure 5:  Distributions in in the $3\ell~{}\PW\PZ$ CR (left) and in
$m_{4\ell}$ in the $4\ell~{}\PZ\PZ$ CR (right). The distributions expected for
the $\PW\PZ$ and $\PZ\PZ$ as well as for other background processes are shown
for the values of nuisance parameters obtained from the ML fit described in
Section 0.9.

The modeling of the reducible $\ell$/misidentification background is validated
in two further CRs, the “$2\ell\text{ss}$ CR” and the “$2\ell+2\tauh$ CR”.
They are based on the SRs of the $2\ell\text{ss}$ and $2\ell+2\tauh$
categories. In the $2\ell\text{ss}$ CR, no $\PQb$ jet veto is applied, and at
least one small-radius jet passing the medium $\PQb$ jet identification is
required. The $2\ell+2\tauh$ CR differs from the SR of the $2\ell+2\tauh$
category in that the sum of $\ell$ plus charges is required to be non-zero,
and no boson veto is applied. The $2\ell\text{ss}$ CR is dominated by events
with misidentified $\ell$, while the $2\ell+2\tauh$ CR is dominated by events
with misidentified . Distributions in the transverse mass in the
$2\ell\text{ss}$ CR and in the mass of the $\PH\PH$ candidate in the
$2\ell+2\tauh$ CR, reconstructed by the algorithm described in Ref. [112], are
shown in Fig. 6. The transverse mass in the $2\ell\text{ss}$ CR is computed
using the leading $\ell$. The data agree well with the background prediction
in both CRs.

Figure 6:  Distributions in in the $2\ell\text{ss}$ CR (left) and in the mass
of the $\PH\PH$ candidate in the $2\ell+2\tauh$ CR (right). The distributions
expected for the misidentified $\ell$/background as well as for other
background processes are shown for the values of nuisance parameters obtained
from the background-only ML fit, in which the $\PH\PH$ signal is constrained
to be zero.

Simulated events are considered as irreducible background if every selected ,
, and candidate matches a prompt MC generator-level counterpart. Events with
at least one selected electron from a photon conversion, and the remaining
$\ell$ and candidates matched to prompt leptons in MC simulation, are
classified as conversion background. Electrons that are misidentified as , and
that are misidentified as are also modeled using the MC simulation. All other
simulated events are discarded, as the $\ell$/misidentification and charge
misidentification backgrounds are estimated from data, as described below.

###  0.7.1 Estimation of the misidentification background

The background from events with misidentified $\ell$ and candidates is
estimated using the “fake factor” or “FF” method from Ref. [110]. An estimate
of this background’s contribution to the SR of each search category is
obtained by selecting a sample of events that satisfy all selection criteria
of the SR for the respective search category, except that the , , and are
required to pass the medium selections instead of the tight ones. The sample
of events thus obtained is referred to as the application region (AR) of the
FF method. Events in which every $\ell$ and satisfies the tight selections are
excluded from the AR.

The prediction for misidentification backgrounds in the SR is obtained by
applying suitably chosen weights $w$ to the events selected in the AR, where
$w$ is given by the expression

$w=(-1)^{n+1}\,\prod_{i=1}^{n}\,\frac{f_{i}(\pt,\eta)}{1-f_{i}(\pt,\eta)}.$
(1)

The product extends over all , , and that pass the medium, but fail the tight
selection criteria, and $n$ refers to the total number of such $\ell$ and .
The symbol $f_{i}(\pt,\eta)$ corresponds to the probability for a single , ,
or that passes the medium selection to also pass the tight one. These
probabilities are measured separately for , , and candidates, and are
parametrized as a function of and $\eta$. The contributions of irreducible
backgrounds to the AR are subtracted based on the MC expectation of such
processes. The alternating sign in Eq. (1) is necessary to avoid double-
counting arising from events with more than one misidentified $\ell$ or [110].

The probabilities $f_{i}(\pt,\eta)$ for electrons and muons are measured in
multijet events, as described in Ref. [98]. The $f_{i}(\pt,\eta)$ for are
measured using $\PZ\to\PGm\PGm$+jets events, where the misidentified
candidates arise from quark or gluon jets. These events are selected by
requiring a muon pair passing the tight selection, with opposite charge and
invariant mass $60<m_{\PGm\PGm}<120$, plus at least one candidate that passes
the medium selection. The leading and subleading muons must have $\pt>25$ and
15, respectively. Events must also pass the $\PQb$ jet veto described in
Section 0.5 to remove background.

###  0.7.2 Estimation of the charge misidentification background

The electron charge misidentification background in the $2\ell\text{ss}$
category is estimated using the method described in Ref. [98]. A sample of
dielectron events passing all selection criteria of the SR of the
$2\ell\text{ss}$ category, except that both electrons are required to have
opposite- instead of same-sign charge, is selected and assigned appropriately
chosen weights. The weights are computed by summing the probabilities for the
charge of either electron to be mismeasured. The probability for the
mismeasurement of the electron charge is determined using $\PZ\to\Pe\Pe$
events, and ranges from under 0.1% in the barrel up to 0.4% in the endcap. The
probability for mismeasuring the charge of muons is negligible [98].

##  0.8 Systematic uncertainties

Multiple sources of systematic uncertainty affect the predicted event yields,
the distributions in the output of the BDT classifiers, or both. These
uncertainties may be theoretical, affecting the predicted cross section or
decay kinematics of the collision process, or experimental, accounting for
differences in object reconstruction and calibration between data and the MC
simulation, as well as for uncertainties on the estimates of the
$\ell$/misidentification and electron charge misidentification background
obtained from data. The systematic uncertainties may be correlated or
uncorrelated across the three data-taking years, and among the various signal
and background processes considered in the analysis.

The SM prediction for the $\Pg\Pg\PH\PH$ production cross section at
$\sqrt{s}=13$has a relative uncertainty of $+6.7$%/$-23.2$%, while the
$\PQq\PQq\PH\PH$ cross section uncertainty is $\pm 2.3\%$. The predicted boson
decay branching fractions to $\PW\PW$, $\PGt\PGt$, and $\PZ\PZ$ have relative
uncertainties of 1.54, 1.65, and 1.54%, respectively [116]. Alternate $\PH\PH$
predictions are generated with the renormalization and factorization scales
varied up and down by a factor of 2. Variations that increase the
factorization scale and decrease the renormalization scale (and vice versa)
are excluded, following the recommendation of Ref. [116]. All theoretical
uncertainties in the $\PH\PH$ signal model are correlated across all three
data-taking years and among the seven search categories. The uncertainties in
the boson decay branching fractions and the effect of renormalization and
factorization scale uncertainties in the signal acceptance impact the
measurement of cross sections for nonresonant as well as resonant $\PH\PH$
production. Conversely, the uncertainties in the SM prediction for the
$\Pg\Pg\PH\PH$ and $\PQq\PQq\PH\PH$ cross sections only affect the measurement
of the $\PH\PH$ production rate as a ratio to the SM prediction.

Theoretical uncertainties also affect the irreducible background prediction.
The relative uncertainties in the cross sections of the dominant $\PW\PZ$ and
$\PZ\PZ$ backgrounds amount to 2.1 and 6.3%, respectively [117, 118, 119]. The
uncertainties in the cross section for the subdominant single boson
backgrounds range from 2 to 9% for $\Pg\Pg\PH\PH$, $\PQq\PQq\PH\PH$, $\PW\PH$,
and $\PZ\PH$. The cross sections for the production of , , or bosons with one
or two top quarks are known with uncertainties of 8–15%. The event yields of
extremely rare backgrounds not mentioned above (, triple boson or four top
quark production) are given a conservative uncertainty of 50%, since the
analysis has little sensitivity to these processes. Following Ref. [98],
background contributions arising from photon conversions are assigned a 30%
yield uncertainty. The theoretical uncertainties affecting background cross
sections are partially correlated among different processes. Here,
contributions arising from uncertainties in the proton PDFs are correlated
among processes with a similar initial state. Processes involving single boson
production are an exception. These uncertainties are uncorrelated from other
background processes but correlated among each other depending on the initial
state. Uncertainties arising from the choice of the renormalization and
factorization scales are correlated for processes with similar final states,
for example among all processes involving diboson production. Uncertainties in
$\alpha_{s}$ are correlated among all background processes. The theoretical
cross section uncertainties for signal processes are uncorrelated with those
of background processes, but otherwise follow the same uncertainty scheme for
proton PDF, scale, and $\alpha_{s}$ contributions. All theoretical cross
section uncertainties are treated as correlated across the different data-
taking years and among all seven search categories.

The rate of the misidentified $\ell$/and electron charge misidentification
backgrounds are assigned a 30% uncertainty in all search categories, to
account for variations in the misidentification rates between the ARs of the
FF method and the multijet ($\PZ\to\PGm\PGm$+jets) event samples used to
measure the $f_{i}(\pt,\eta)$ for and (). In the $3\ell+1\tauh$ and
$1\ell+3\tauh$ categories, an additional uncertainty of 30% (uncorrelated with
the other 30% uncertainty) is assigned to the rate of the misidentified
$\ell$/background, to account for the extra uncertainty arising from the
modified selection criteria that suppress the misidentification of electrons
as . The effect of statistical uncertainties in the probabilities
$f_{i}(\pt,\eta)$ for electrons and muons is evaluated by varying these
probabilities in bins of and $\eta$ and determining the resulting change in
the shape of the BDT classifier output distribution obtained for the
misidentified $\ell$/background. For , the effect of statistical uncertainties
in $f_{i}(\pt,\eta)$ is evaluated by fitting the probabilities in bins of
$\eta$ with functions that are linear in , varying the slope of these
functions up and down within the uncertainties obtained from the fit, and
determining the resulting change in the shape of the BDT classifier output
distribution.

An additional uncertainty in the BDT output shape in each category is
evaluated for events with a nonprompt or misidentified $\ell$ or as follows:
Simulated events passing all signal selection criteria are compared to those
with at least one $\ell$ or candidate failing the tight identification
criteria, scaled according to the FF method described in Section 0.7, but with
the probabilities $f_{i}(\pt,\eta)$ taken from the MC simulation instead of
from the data. The ratio of these two shapes is fitted with a linear function,
which is convoluted with the misidentified $\ell$/background prediction from
the data to serve as an uncertainty in the BDT output shape for these events
in the SR. The systematic uncertainties associated with the misidentified
$\ell$/background prediction and the uncertainty associated with the electron
charge misidentification rate are treated as uncorrelated among the different
data-taking years.

Uncertainties in the modeling of the trigger and object reconstruction
efficiency affect all signal and background processes that are estimated using
MC simulation. Trigger efficiencies for events with at least two $\ell$ are
compared between data and MC simulation in control regions enriched in the ,
$\PW\PZ$, and $\PZ\PZ$ background processes, as a function of lepton flavor, ,
and $\eta$. This results in a small -dependent uncertainty correlated between
the $2\ell\text{ss}$ and $2\ell+2\tauh$ categories, and a 1% normalization
uncertainty, which is correlated among the $3\ell$, $3\ell+1\tauh$, and
$4\ell$ categories. The data-to-simulation agreement in trigger efficiency for
the $4\tauh$ and $1\ell+3\tauh$ categories is computed using an independent
set of data, as a function of the and $\eta$ of the $\ell$ and all , and the
reconstructed decay modes of all . The trigger uncertainties for these two
categories are treated as uncorrelated. All systematic uncertainties related
to trigger modeling are correlated across different physics processes, but
uncorrelated among the three data-taking years.

The uncertainties in the reconstruction and identification efficiencies for ,
, and candidates have been measured in boson enriched regions in data for each
level of identification criteria (tight, medium, and loose), and are applied
to each event as a function of and $\eta$ for leptons and of and the
reconstructed hadronic decay mode for . The reconstructed energy has an
uncertainty of around 1%, depending on the data-taking year and reconstructed
decay mode. These uncertainties affect the predicted rate and BDT output shape
for signal and background, and are correlated among the different physics
processes, but uncorrelated across different data-taking years.

The jet energy scale and resolution are determined using dijet control regions
[106, 120]. The jet energy scale is evaluated using 11 separate components,
accounting for partial correlations between the data recorded in different
years. The jet energy resolution uncertainty is uncorrelated among the three
data-taking years. Jet energy uncertainties are also propagated to the
calculation. An additional uncertainty in comes from uncertainty in the energy
of “unclustered” PF hadrons (PF hadrons not clustered into either small- or
large-radius jets), which is uncorrelated across different years. The
probability for true $\PQb$ jets to fail the multivariate $\PQb$ jet
identification criteria, or for jets from gluons or light flavored quarks to
be misidentified as $\PQb$ jets, is compared in data and MC simulation in
event regions that are enriched in light-flavor quark or gluon, or heavy-
flavor jets. The resulting uncertainty in the data-to-simulation agreement
affects the yields and BDT output shapes of multiple physics processes. The
statistical component of this uncertainty is treated as uncorrelated across
different data-taking years, while other experimental sources are correlated.

The integrated luminosities for data collected in 2016, 2017, and 2018 have
1.2–2.5% individual uncertainties [55, 56, 57], while the overall uncertainty
for the 2016–2018 period is 1.6%. The uncertainty in the measured cross
section for inelastic $\Pp\Pp$ collisions, amounting to 5% [121], is taken
into account by varying the number of pileup interactions in MC simulation,
which impacts the jet reconstruction and the isolation of $\ell$ and .

The sources of systematic uncertainty which create the largest uncertainties
in the measured ratio of the $\PH\PH$ production cross section to its SM
prediction are the theoretical uncertainties in the $\PH\PH$ production cross
section and decay branching fractions (25%), the uncertainties in the rate and
shape of backgrounds from misidentified $\ell$ or (22%), and in the rates of
backgrounds modeled using MC simulation (13%). These uncertainties in the
signal measurement are determined by removing uncertainties that correspond to
a given systematic source from an ML fit to pseudodata, as described in
Section 0.9, and subtracting the obtained uncertainty in the signal
measurement in quadrature from the total uncertainty. The impacts of
systematic uncertainties are small compared to the effect of the statistical
uncertainty in the data (79%), and are comparable to the statistical
uncertainties in the distributions in the BDT classifier output for background
processes (33%). The latter includes the effect of statistical uncertainties
in the MC simulation as well as in the $\ell$/misidentification and electron
charge misidentification backgrounds obtained from data. All other sources of
uncertainty have an impact of 5% or less.

##  0.9 Results

The data selected in the seven search categories are tested against multiple
$\PH\PH$ production hypotheses: the SM prediction; variations of the SM
coupling strength modifiers $\kappa_{\lambda}$, $\kappa_{\PQt}$,
$\kappa_{\PV}$, and $\kappa_{2\PV}$; the effective couplings $\text{c}_{\Pg}$,
$\text{c}_{2\Pg}$, and $\text{c}_{2}$ in the EFT approach; and resonant
production of boson pairs originating from the decay of heavy particles with
spins of 0 or 2 and masses $m_{\textrm{X}}$ ranging from 250 to 1000. In each
case, the entire data set is fit simultaneously to a model composed of the
background prediction (with uncertainties) and the $\PH\PH$ signal hypothesis
under consideration. The distributions in in the $3\ell~{}\PW\PZ$ CR and in
$m_{4\ell}$ in the $4\ell~{}\PZ\PZ$ CR, enriched in $\PW\PZ$ and $\PZ\PZ$
backgrounds and shown in Fig. 5, are included in each of these fits, in order
to obtain in-situ constraints on the systematic uncertainties described in
Section 0.8. This in turn reduces the uncertainties in signal and background
predictions in the seven search categories.

The SM “signal strength” parameter $\mu$ is defined as the ratio of the
measured $\PH\PH$ production cross section to its predicted value in the SM.
This parameter modifies the expected signal yield by the same proportion in
each category. By contrast, variations in the $\kappa$ modifiers may affect
the signal yields in each category differently, and also change the BDT
classifier output shape for $\PH\PH$ events. The twenty benchmark scenarios
spanning combinations of $\kappa_{\lambda}$, $\kappa_{\PQt}$,
$\text{c}_{\Pg}$, $\text{c}_{2\Pg}$, and $\text{c}_{2}$ values in the coupling
parameter space each correspond to different kinematic distributions, so the
$\PH\PH$ production cross section for each point is measured separately.
Similarly, signal efficiency and BDT classifier output shapes vary
dramatically for different resonant masses, and thus a separate measurement is
performed for each mass and spin hypothesis. The SM signal strength
measurement is performed using the output of the BDT classifier that has been
trained for SM nonresonant $\PH\PH$ production, while the $\kappa_{\lambda}$
measurement uses the BDT trained for benchmark scenario JHEP04 BM7. In the
scenario JHEP04 BM7, the $m_{\PH\PH}$ value tends to be close to the lower
limit of 250, which matches the event kinematics for non-resonant $\PH\PH$
production in the $\kappa_{\lambda}$ range in which we expect to set the
limit. When setting limits on the twenty different benchmark scenarios, the
binary BDT inputs correspond to the given scenario, or in case of the
benchmarks from Ref. [58] the kinematically closest scenario. In case of
resonant $\PH\PH$ production, the BDT input for the resonance mass is set to
the $m_{\textrm{X}}$ value for which the limit is computed.

The SM signal strength is measured using a profile likelihood test statistic
[122], with systematic uncertainties treated as nuisance parameters $\theta$
in a frequentist approach [123]. The effect of variations in $\theta$ on the
shape of the BDT classifier output distribution for the $\PH\PH$ signal and
for background processes is incorporated into the ML fit using the technique
described in Ref. [124]. Statistical uncertainties in these distributions are
also taken into account using the approach detailed in Ref. [124]. The
likelihood ratio $q_{\mu}$ for a fixed “test” signal strength value $\mu$ is

$\displaystyle q_{\mu}$
$\displaystyle=-2\Delta\ln\mathcal{L}=-2\ln\frac{\mathcal{L}(\text{data}|\mu,\hat{\theta}_{\mu})}{\mathcal{L}(\text{data}|\hat{\mu},\hat{\theta})},$

where $\hat{\mu}$ and $\hat{\theta}$ are the signal strength and nuisance
parameter values that give the maximum value of the likelihood function
$\mathcal{L}$ for the given set of data (requiring $\hat{\mu}\geq 0$), and
$\hat{\theta}_{\mu}$ is the set of $\theta$ values which maximize
$\mathcal{L}$ for the fixed $\mu$. The 95% confidence level () upper limit for
$\mu$ is obtained using the criterion [125, 126], with $q_{\mu}$ set to $0$
when $\mu<\hat{\mu}$. The probabilities to observe a given value of the
likelihood ratio $q_{\mu}$ under the signal-plus-background and background-
only hypotheses are computed using the asymptotic approximation from Ref.
[122]. The limits on $\mu$ obtained using the asymptotic approximation, with
all seven search categories and the two CRs for the $\PW\PZ$ and $\PZ\PZ$
backgrounds fit simultaneously, match the limits obtained with toy MC
experiments [123] within 10%. The SM coupling strength modifiers and the cross
sections for the various $\PH\PH$ production hypotheses are measured by
scanning the likelihood ratio $q_{\mu}$ as a function of $\mu$. Theoretical
and experimental uncertainties affecting the signal and background yields or
the shape of the BDT classifier output distributions may be correlated or
uncorrelated across different years, search categories, and BDT output bins,
as described in Section 0.8.

The fits to all BDT output distributions in the data are consistent with the
SM background prediction in all search categories, within statistical and
systematic uncertainties. Distributions in the output of the BDT classifier
for SM nonresonant $\PH\PH$ production in the seven search categories are
shown in Figs. 7 and 8. The data excess in the rightmost bin of the BDT
classifier output distribution for the $3\ell$ category is not statistically
significant: 11 events are observed in this bin, while $5.2\pm 0.7$ are
expected from background processes, amounting to a local significance of about
1.7 standard deviations.

For the case of nonresonant $\PH\PH$ production with event kinematics as
predicted by the SM, the best-fit value of the $\PH\PH$ production rate,
obtained from the simultaneous fit of all seven search categories with the
$3\ell~{}\PW\PZ$ and $4\ell~{}\PZ\PZ$ CRs, amounts to $\hat{\mu}=2\pm 10$
times the SM expectation. The measured value of the signal strength refers to
the sum of $\Pg\Pg\PH\PH$ and $\PQq\PQq\PH\PH$ production and is compatible
with the SM as well as with the background-only hypothesis. The observed
(expected) 95% upper limit on the cross section for nonresonant $\PH\PH$
production is 651 (592) fb. Taking into account the theoretical uncertainties
in the SM $\PH\PH$ production cross section, this corresponds to an observed
(expected) limit on the nonresonant $\PH\PH$ production rate of 21.3 (19.4)
times the SM expectation. These limits are shown in Fig. 9, for individual
categories and for the result obtained from the combination of all seven
search categories. The $3\ell~{}\PW\PZ$ and $4\ell~{}\PZ\PZ$ CRs are included
in all of these fits. We refer to the combined result as the
“$\PH\PH\to\text{multilepton}$” result. The $3\ell$ and $1\ell+3\tauh$
categories are the most sensitive to SM $\PH\PH$ production, followed closely
by the other categories.

Figure 7:  Distribution in the output of the BDT trained for nonresonant
$\PH\PH$ production and evaluated for the benchmark scenario JHEP04 BM7 for
the $2\ell\text{ss}$ (upper left), $3\ell$ (upper right), and $4\ell$ (lower)
categories. The SM $\PH\PH$ signal is shown for a cross section amounting to
30 times the value predicted in the SM. The distributions expected for the
background processes are shown for the values of nuisance parameters obtained
from the ML fit of the signal+background hypothesis to the data.

Figure 8:  Distribution in the output of the BDT trained for nonresonant
$\PH\PH$ production and evaluated for the benchmark scenario JHEP04 BM7 for
the $3\ell+1\tauh$ (upper left), $2\ell+2\tauh$ (upper right), $1\ell+3\tauh$
(lower left), and $4\tauh$ (lower right) categories. The SM $\PH\PH$ signal is
shown for a cross section amounting to 30 times the value predicted in the SM.
The distributions expected for the background processes are shown for the
values of nuisance parameters obtained from the ML fit of the
signal+background hypothesis to the data.  Figure 9:  Observed and expected
95% upper limits on the SM $\PH\PH$ production cross section, obtained for
both individual search categories and from a simultaneous fit of all seven
categories combined.

The observed (expected) 95% interval for the boson trilinear self-coupling
strength modifier is measured to be $-6.9<\kappa_{\lambda}<11.1$
($-6.9<\kappa_{\lambda}<11.7$). The upper limit on $\kappa_{\lambda}$ is one
of the strongest constraints on this fundamental SM parameter to date, with
only $\PH\PH$ searches in the $\PQb\PQb\PGg\PGg$ [41, 42] and
$\PQb\PQb\PQb\PQb$ [43] decay modes providing greater sensitivity. The
observed and expected upper limits on the $\PH\PH$ production cross section as
a function of $\kappa_{\lambda}$, obtained from the simultaneous fit of all
seven search categories, are shown in Fig. 10, along with the limits obtained
for each category individually. The $3\ell~{}\PW\PZ$ and $4\ell~{}\PZ\PZ$ CRs
are again included in all of these fits.

Figure 10:  Observed and expected 95% upper limits on the $\PH\PH$ production
cross section as a function of the boson self-coupling strength modifier
$\kappa_{\lambda}$. All boson couplings other than $\lambda$ are assumed to
have the values predicted in the SM. The left plot shows the result obtained
by combining all seven search categories, while the right plot shows the
limits obtained for each category separately. The red curve in the left plot
represents the SM prediction for the $\PH\PH$ production cross section as a
function of $\kappa_{\lambda}$, and the red shaded band the theoretical
uncertainty in this prediction.

The observed and expected limits on the $\Pg\Pg\PH\PH$ production cross
section for the twenty benchmark scenarios are shown in Fig. 11 and summarized
in Table 0.9. The $\PQq\PQq\PH\PH$ process can be safely neglected for these
measurements. The observed (expected) limits on nonresonant $\PH\PH$
production in the different benchmark scenarios range from 0.21 to 1.09 (0.16
to 1.16) pb, depending on the scenario. These limits are comparable to, but
somewhat higher, than those obtained by the CMS measurement in the
$\PQb\PQb\PGg\PGg$ final state [41]. The variation in expected limits reflects
differences in the $m_{\PH\PH}$ distribution among the benchmark scenarios,
which in turn affect the and angles between the particles produced in the
boson decays. As a consequence, the signal acceptance can change, as well as
the separation of the $\PH\PH$ signal from backgrounds through the BDT
classifiers described in Section 0.6. The most and least stringent limits on
the cross section are expected for the benchmark scenarios JHEP04 BM2 and BM7,
respectively. The former has a pronounced tail of the $m_{\PH\PH}$
distribution extending to high values, while the latter is characterized by
low $m_{\PH\PH}$ values, as seen in Fig. 5 of Ref. [24].

Figure 11:  Observed and expected 95% upper limits on the $\PH\PH$ production
cross section for the twelve benchmark scenarios from Ref. [24], the
additional benchmark scenario 8a from Ref. [59], the seven benchmark scenarios
from Ref. [58], and for the SM. The upper plot shows the result obtained by
combining all seven search categories, while the lower plot shows the limits
obtained for each category separately, and the combined limit.

Observed (expected) 95% upper limits on the $\Pg\Pg\PH\PH$ production cross
section for the twelve benchmark scenarios from Ref. [24], the additional
benchmark scenario 8a from Ref. [59] and the seven benchmark scenarios from
Ref. [58]. The corresponding observed (expected) upper limit for the SM is 652
(583) fb. The limits correspond to the combination of all seven search
categories.  JHEP04 Observed (expected) benchmark limit [fb] BM1 469 (354) BM2
205 (159) BM3 563 (447) BM4 677 (600) BM5 439 (263) BM6 739 (584) BM7 1090
(1156) BM8 495 (336) BM9 541 (298) BM10 988 (855) BM11 795 (572) BM12 897
(898) BM8a 608 (353) JHEP03 Observed (expected) benchmark limit [fb] BM1 888
(650) BM2 828 (632) BM3 538 (293) BM4 559 (436) BM5 556 (313) BM6 660 (518)
BM7 525 (280)

Figure 12 shows the observed and expected upper limits on the $\PH\PH$
production cross section as a function of the coupling $\text{c}_{2}$, and the
region excluded in the $\kappa_{\PQt}$–$\text{c}_{2}$ plane. The effects of
variations in $\kappa_{\lambda}$ and $\kappa_{\PQt}$ on the rate of the SM
single boson background [21] and on the boson decay branching fractions [20]
are taken into account when computing these limits and those shown in Fig. 10.
The magnitude of these effects is typically 5 to 10% within the scanned range
of $\kappa_{\lambda}$ and $\kappa_{\PQt}$. Assuming $\kappa_{\PQt}$ and
$\kappa_{\lambda}$ are both equal to 1, the coupling $\text{c}_{2}$ is
observed (expected) to be constrained to the interval
$-1.05<\text{c}_{2}<1.48$ ($-0.96<\text{c}_{2}<1.37$) at 95% .

Figure 12:  Observed and expected limits on the $\PH\PH$ production cross
section as a function of the effective coupling $\text{c}_{2}$ (left), and the
region excluded in the $\kappa_{\PQt}$–$\text{c}_{2}$ plane (right). All
limits are computed at 95% . boson couplings other than the ones shown in the
plots ($\text{c}_{2}$ in the left plot and $\text{c}_{2}$ and $\kappa_{\PQt}$
in the right plot) are assumed to have the values predicted by the SM.

Figure 13 shows the observed and expected limits on the resonant $\PH\PH$
production cross section as a function of $m_{\textrm{X}}$ for a spin-0 or
spin-2 particle X decaying to $\PH\PH$. The mass points probed are listed in
the fourth paragraph of Section 0.3. The limits are expected to become more
stringent as $m_{\textrm{X}}$ increases, as the acceptance for the $\PH\PH$
signal increases and the signal can be more easily distinguished from
backgrounds. The observed (expected) 95% upper limits on the resonant $\PH\PH$
production cross section range from 0.18 to 0.90 (0.08 to 1.06) pb, depending
on the mass and spin. Compared to existing combinations of CMS and ATLAS
results [50, 51], this analysis has superior sensitivity at very low masses
(close to 250), owing to the efficient reconstruction and identification of
low-leptons in the CMS detector. The recent ATLAS search in the
$\PQb\PQb\PGg\PGg$ final state achieves somewhat lower limits all the way down
to 250 [42]. For $m_{\textrm{X}}\gtrsim 600$, the observed limit is less
stringent than the expected limit, due to a small excess of events in the data
that is concentrated near $m_{\textrm{X}}=750$in the $2\ell\text{ss}$ and
$3\ell$ categories. The distributions in the output of the BDT classifier
targeting resonances with spin 2 and mass 750in the $2\ell\text{ss}$ and
$3\ell$ categories are shown in Fig. 14. A small excess of events can be seen
in the rightmost bin of both distributions. In the $2\ell\text{ss}$ ($3\ell$)
category, 42 (17) events are observed in this bin in the data, while $27.3\pm
2.9$ ($8.0\pm 1.0$) are expected from background processes, amounting to a
local significance of about 2.1 (2.1) standard deviations. The excess affects
the observed limits in a broad mass range from 600 to 1000. The presence of
multiple neutrinos in $\PH\PH$ signal events in these categories, coming from
boson or lepton decays, limits the experimental resolution on $m_{\textrm{X}}$
and causes the BDT classifier output distributions to be highly correlated for
resonances of similar mass. No significant excess is observed in any of the
other five search categories. Tabulated results are provided in the HEPData
record for this analysis [127].

Figure 13:  Observed and expected 95% upper limits on the production of new
particles X of spin 0 (upper) and spin 2 (lower) and mass $m_{\textrm{X}}$ in
the range 250–1000, which decay to boson pairs. The plot on the left shows the
result obtained by combining all seven search categories, while the plot on
the right shows the limits obtained for each category separately, and the
combined limit.

Figure 14:  Distribution in BDT classifier output for resonances of spin 2 and
mass 750in the $2\ell\text{ss}$ (left) and $3\ell$ (right) categories. The
resonant $\PH\PH$ signal is shown for a cross section amounting to 1 pb. The
distributions expected for the background processes are shown for the values
of nuisance parameters obtained from the ML fit of the signal+background
hypothesis to the data.

##  0.10 Summary

The results of a search for nonresonant as well as resonant Higgs boson pair
($\PH\PH$) production in final states with multiple reconstructed leptons,
including electrons and muons ($\ell$) as well as hadronically decaying tau
leptons (), has been presented. The search targets the $\PH\PH$ decay modes
$\PW\PW\PW\PW$, $\PW\PW\Pgt\Pgt$, and $\Pgt\Pgt\Pgt\Pgt$, using proton-proton
collision data recorded by the CMS experiment at a center-of-mass energy of
13and corresponding to an integrated luminosity of 138. Seven search
categories, distinguished by $\ell$ and multiplicity, are included in the
analysis: $2\ell\text{ss}$, $3\ell$, $4\ell$, $3\ell+1\tauh$, $2\ell+2\tauh$,
$1\ell+3\tauh$, and $4\tauh$, where “ss” indicates an $\ell\ell$ pair with the
same charge. No evidence for a signal is found in the data. Upper limits on
the cross section for nonresonant as well as resonant $\PH\PH$ production are
set. The observed (expected) limits on the nonresonant $\PH\PH$ production
cross section in twenty EFT benchmark scenarios range from 0.21 to 1.09 (0.16
to 1.16) pb at 95% confidence level (), depending on the scenario. For
nonresonant $\PH\PH$ production with event kinematics as predicted by the
standard model (SM), the observed (expected) 95% upper limit on the $\PH\PH$
production rate is 21.3 (19.4) times the rate expected in the SM. The results
of the search for nonresonant $\PH\PH$ production are used to exclude regions
in the plane of the boson coupling to the top quark, $\text{y}_{\PQt}$, and of
the trilinear Higgs boson self-coupling, $\lambda$. Assuming $\text{y}_{\PQt}$
has the value expected in the SM, the observed (expected) 95% interval for
$\lambda$ is between $-6.9$ and 11.1 ($-6.9$ and 11.7) times the value
expected in the SM. The resonant production of boson pairs, resulting from
decays of new heavy particles X with mass $m_{\textrm{X}}$, is probed within
the mass range 250–1000. The corresponding observed (expected) 95% upper
limits on the cross section for resonant $\PH\PH$ production range from 0.18
to 0.90 (0.08 to 1.06) pb, depending on the mass and spin of the resonance.

## References

  * [1] ATLAS Collaboration, “Observation of a new particle in the search for the standard model Higgs boson with the ATLAS detector at the LHC”,  Phys. Lett. B 716 (2012) 1, 10.1016/j.physletb.2012.08.020, arXiv:1207.7214.
  * [2] CMS Collaboration, “Observation of a new boson at a mass of 125 with the CMS experiment at the LHC”,  Phys. Lett. B 716 (2012) 30, 10.1016/j.physletb.2012.08.021, arXiv:1207.7235.
  * [3] CMS Collaboration, “Observation of a new boson with mass near 125 in $\Pp\Pp$ collisions at $\sqrt{s}=7$ and 8”,  JHEP 06 (2013) 081, 10.1007/JHEP06(2013)081, arXiv:1303.4571.
  * [4] CMS Collaboration, “Combined measurements of Higgs boson couplings in proton-proton collisions at $\sqrt{s}=13\TeV$”,  Eur. Phys. J. C 79 (2019) 421, 10.1140/epjc/s10052-019-6909-y, arXiv:1809.10733.
  * [5] ATLAS Collaboration, “Combined measurements of Higgs boson production and decay using up to $80\fbinv$ of proton-proton collision data at $\sqrt{s}=13\TeV$ collected with the ATLAS experiment”,  Phys. Rev. D 101 (2020) 012002, 10.1103/PhysRevD.101.012002, arXiv:1909.02845.
  * [6] CMS Collaboration, “First evidence for off-shell production of the Higgs boson and measurement of its width”, 2022. arXiv:2202.06923. Submitted to Nature Physics.
  * [7] F. Englert and R. Brout, “Broken symmetry and the mass of gauge vector mesons”,  Phys. Rev. Lett. 13 (1964) 321, 10.1103/PhysRevLett.13.321.
  * [8] P. W. Higgs, “Broken symmetries, massless particles and gauge fields”,  Phys. Lett. 12 (1964) 132, 10.1016/0031-9163(64)91136-9.
  * [9] P. W. Higgs, “Broken symmetries and the masses of gauge bosons”,  Phys. Rev. Lett. 13 (1964) 508, 10.1103/PhysRevLett.13.508.
  * [10] G. S. Guralnik, C. R. Hagen, and T. W. B. Kibble, “Global conservation laws and massless particles”,  Phys. Rev. Lett. 13 (1964) 585, 10.1103/PhysRevLett.13.585.
  * [11] P. W. Higgs, “Spontaneous symmetry breakdown without massless bosons”,  Phys. Rev. 145 (1966) 1156, 10.1103/PhysRev.145.1156.
  * [12] T. W. B. Kibble, “Symmetry breaking in non-Abelian gauge theories”,  Phys. Rev. 155 (1967) 1554, 10.1103/PhysRev.155.1554.
  * [13] S. L. Glashow, “Partial symmetries of weak interactions”,  Nucl. Phys. 22 (1961) 579, 10.1016/0029-5582(61)90469-2.
  * [14] S. Weinberg, “A model of leptons”,  Phys. Rev. Lett. 19 (1967) 1264, 10.1103/PhysRevLett.19.1264.
  * [15] A. Salam, “Weak and electromagnetic interactions”,  Conf. Proc. C 680519 (1968) 367, 10.1142/9789812795915_0034.
  * [16] D. de Florian, I. Fabre, and J. Mazzitelli, “Triple Higgs production at hadron colliders at NNLO in QCD”,  JHEP 03 (2020) 155, 10.1007/JHEP03(2020)155, arXiv:1912.02760.
  * [17] O. Aberle et al., “High-Luminosity Large Hadron Collider (HL-LHC): Technical design report V. 0.1”,  CERN Yellow Rep. Monogr. 4 (2017) 10.23731/CYRM-2017-004.
  * [18] M. Grazzini et al., “Higgs boson pair production at NNLO with top quark mass effects”,  JHEP 05 (2018) 059, 10.1007/JHEP05(2018)059, arXiv:1803.02463.
  * [19] F. A. Dreyer and A. Karlberg, “Vector-boson fusion Higgs pair production at N3LO”,  Phys. Rev. D 98 (2018) 114016, 10.1103/PhysRevD.98.114016, arXiv:1811.07906.
  * [20] G. Degrassi, P. P. Giardino, F. Maltoni, and D. Pagani, “Probing the Higgs self coupling via single Higgs production at the LHC”,  JHEP 12 (2016) 080, 10.1007/JHEP12(2016)080, arXiv:1607.04251.
  * [21] F. Maltoni, D. Pagani, A. Shivaji, and X. Zhao, “Trilinear Higgs coupling determination via single-Higgs differential measurements at the LHC”,  Eur. Phys. J. C 77 (2017) 887, 10.1140/epjc/s10052-017-5410-8, arXiv:1709.08649.
  * [22] W. Buchmuller and D. Wyler, “Effective Lagrangian analysis of new interactions and flavor conservation”,  Nucl. Phys. B 268 (1986) 621, 10.1016/0550-3213(86)90262-2.
  * [23] B. Grzadkowski, M. Iskrzynski, M. Misiak, and J. Rosiek, “Dimension-six terms in the standard model Lagrangian”,  JHEP 10 (2010) 085, 10.1007/JHEP10(2010)085, arXiv:1008.4884.
  * [24] A. Carvalho et al., “Higgs pair production: choosing benchmarks with cluster analysis”,  JHEP 04 (2016) 126, 10.1007/JHEP04(2016)126, arXiv:1507.02245.
  * [25] N. Craig, J. Galloway, and S. Thomas, “Searching for signs of the second Higgs doublet”, 2013. arXiv:1305.2424.
  * [26] D. T. Nhung, M. Mühlleitner, J. Streicher, and K. Walz, “Higher order corrections to the trilinear Higgs self-couplings in the real NMSSM”,  JHEP 11 (2013) 181, 10.1007/JHEP11(2013)181, arXiv:1306.3926.
  * [27] R. Grober and M. Mühlleitner, “Composite Higgs boson pair production at the LHC”,  JHEP 06 (2011) 020, 10.1007/JHEP06(2011)020, arXiv:1012.1562.
  * [28] R. Contino et al., “Strong double Higgs production at the LHC”,  JHEP 05 (2010) 089, 10.1007/JHEP05(2010)089, arXiv:1002.1011.
  * [29] C. Englert, T. Plehn, D. Zerwas, and P. M. Zerwas, “Exploring the Higgs portal”,  Phys. Lett. B 703 (2011) 298, 10.1016/j.physletb.2011.08.002, arXiv:1106.3097.
  * [30] J. M. No and M. Ramsey-Musolf, “Probing the Higgs portal at the LHC through resonant di-Higgs production”,  Phys. Rev. D 89 (2014) 095031, 10.1103/PhysRevD.89.095031, arXiv:1310.6035.
  * [31] L. Randall and R. Sundrum, “A large mass hierarchy from a small extra dimension”,  Phys. Rev. Lett. 83 (1999) 3370, 10.1103/PhysRevLett.83.3370, arXiv:hep-ph/9905221.
  * [32] K. Cheung, “Phenomenology of radion in Randall-Sundrum scenario”,  Phys. Rev. D 63 (2001) 056007, 10.1103/PhysRevD.63.056007, arXiv:hep-ph/0009232.
  * [33] CMS Collaboration, “Search for heavy resonances decaying to a pair of Lorentz-boosted Higgs bosons in final states with leptons and a bottom quark pair at $\sqrt{s}$ = 13 TeV”, 2021. arXiv:2112.03161. Submitted to JHEP.
  * [34] ATLAS Collaboration, “Search for pair production of Higgs bosons in the $\PQb\PAQb\PQb\PAQb$ final state using proton-proton collisions at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  JHEP 01 (2019) 030, 10.1007/JHEP01(2019)030, arXiv:1804.06174.
  * [35] U. Baur, T. Plehn, and D. L. Rainwater, “Measuring the Higgs boson self-coupling at the LHC and finite top mass matrix elements”,  Phys. Rev. Lett. 89 (2002) 151801, 10.1103/PhysRevLett.89.151801, arXiv:hep-ph/0206024.
  * [36] U. Baur, T. Plehn, and D. L. Rainwater, “Determining the Higgs boson self-coupling at hadron colliders”,  Phys. Rev. D 67 (2003) 033003, 10.1103/PhysRevD.67.033003, arXiv:hep-ph/0211224.
  * [37] Q. Li, Z. Li, Q.-S. Yan, and X. Zhao, “Probe Higgs boson pair production via the $3\ell\,2\text{j}+\ptmiss$ mode”,  Phys. Rev. D 92 (2015) 014015, 10.1103/PhysRevD.92.014015, arXiv:1503.07611.
  * [38] A. Adhikary et al., “Revisiting the non-resonant Higgs pair production at the HL-LHC”,  JHEP 07 (2018) 116, 10.1007/JHEP07(2018)116, arXiv:1712.05346.
  * [39] J. Ren et al., “LHC search of new Higgs boson via resonant di-Higgs production with decays into $4\PW$”,  JHEP 06 (2018) 090, 10.1007/JHEP06(2018)090, arXiv:1706.05980.
  * [40] ATLAS Collaboration, “Search for Higgs boson pair production in the $\PW\PW^{\ast}\PW\PW^{\ast}$ decay channel using ATLAS data recorded at $\sqrt{s}=13\TeV$”,  JHEP 05 (2019) 124, 10.1007/JHEP05(2019)124, arXiv:1811.11028.
  * [41] CMS Collaboration, “Search for nonresonant Higgs boson pair production in final states with two bottom quarks and two photons in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JHEP 03 (2021) 257, 10.1007/JHEP03(2021)257, arXiv:2011.12373.
  * [42] ATLAS Collaboration, “Search for Higgs boson pair production in the two bottom quarks plus two photons final state in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$ with the ATLAS detector”, 2021. arXiv:2112.11876. Submitted to Phys. Rev. D.
  * [43] CMS Collaboration, “Search for Higgs boson pair production in the four b quark final state in proton-proton collisions at $\sqrt{s}$ = 13 TeV”, 2022. arXiv:2202.09617. Submitted to Phys. Rev. Lett.
  * [44] CMS Collaboration, “Search for resonant pair production of Higgs bosons decaying to bottom quark-antiquark pairs in proton-proton collisions at $13\TeV$”,  JHEP 08 (2018) 152, 10.1007/JHEP08(2018)152, arXiv:1806.03548.
  * [45] CMS Collaboration, “Search for Higgs boson pair production in events with two bottom quarks and two tau leptons in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$”,  Phys. Lett. B 778 (2018) 101, 10.1016/j.physletb.2018.01.001, arXiv:1707.02909.
  * [46] ATLAS Collaboration, “Search for resonant and non-resonant Higgs boson pair production in the $\PQb\PAQb\PGt^{+}\PGt^{-}$ decay channel in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  Phys. Rev. Lett. 121 (2018), no. 19, 191801, 10.1103/PhysRevLett.121.191801, arXiv:1808.00336. [Erratum: 10.1103/PhysRevLett.122.089901].
  * [47] CMS Collaboration, “Search for resonant and non-resonant Higgs boson pair production in the $\PQb\PAQb\ell\PGn\ell\PGn$ final state in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JHEP 01 (2018) 054, 10.1007/JHEP01(2018)054, arXiv:1708.04188.
  * [48] ATLAS Collaboration, “Search for non-resonant Higgs boson pair production in the $\PQb\PAQb\ell\PGn\ell\PGn$ final state with the ATLAS detector in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$”,  Phys. Lett. B 801 (2020) 135145, 10.1016/j.physletb.2019.135145, arXiv:1908.06765.
  * [49] ATLAS Collaboration, “Search for Higgs boson pair production in the $\PGg\PGg\PW\PW^{\ast}$ channel using $\Pp\Pp$ collision data recorded at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  Eur. Phys. J. C 78 (2018) 1007, 10.1140/epjc/s10052-018-6457-x, arXiv:1807.08567.
  * [50] CMS Collaboration, “Combination of searches for Higgs boson pair production in proton-proton collisions at $\sqrt{s}=13\TeV$”,  Phys. Rev. Lett. 122 (2019) 121803, 10.1103/PhysRevLett.122.121803, arXiv:1811.09689.
  * [51] ATLAS Collaboration, “Combination of searches for Higgs boson pairs in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$ with the ATLAS detector”,  Phys. Lett. B 800 (2020) 135103, 10.1016/j.physletb.2019.135103, arXiv:1906.02025.
  * [52] CMS Collaboration, “Performance of the CMS level-1 trigger in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JINST 15 (2020) P10017, 10.1088/1748-0221/15/10/P10017, arXiv:2006.10165.
  * [53] CMS Collaboration, “The CMS trigger system”,  JINST 12 (2017) P01020, 10.1088/1748-0221/12/01/P01020, arXiv:1609.02366.
  * [54] CMS Collaboration, “The CMS experiment at the CERN LHC”,  JINST 3 (2008) S08004, 10.1088/1748-0221/3/08/S08004.
  * [55] CMS Collaboration, “Precision luminosity measurement in proton-proton collisions at $\sqrt{s}=13\TeV$ in 2015 and 2016 at CMS”,  Eur. Phys. J. C 81 (2021) 800, 10.1140/epjc/s10052-021-09538-2, arXiv:2104.01927.
  * [56] CMS Collaboration, “CMS luminosity measurement for the 2017 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-17-004, 2018.
  * [57] CMS Collaboration, “CMS luminosity measurement for the 2018 data-taking period at $\sqrt{s}=13\TeV$”, CMS Physics Analysis Summary CMS-PAS-LUM-18-002, 2019.
  * [58] M. Capozi and G. Heinrich, “Exploring anomalous couplings in Higgs boson pair production through shape analysis”,  JHEP 03 (2020) 091, 10.1007/JHEP03(2020)091, arXiv:1908.08923.
  * [59] G. Buchalla et al., “Higgs boson pair production in non-linear effective field theory with full $m_{\PQt}$-dependence at NLO QCD”,  JHEP 09 (2018) 057, 10.1007/JHEP09(2018)057, arXiv:1806.05162.
  * [60] A. Carvalho et al., “On the reinterpretation of non-resonant searches for Higgs boson pairs”,  JHEP 02 (2021) 049, 10.1007/JHEP02(2021)049, arXiv:1710.08261.
  * [61] J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”,  JHEP 07 (2014) 079, 10.1007/JHEP07(2014)079, arXiv:1405.0301.
  * [62] A. Kalogeropoulos and J. Alwall, “The SysCalc code: A tool to derive theoretical systematic uncertainties”, 1, 2018. arXiv:1801.08401.
  * [63] P. Nason, “A new method for combining NLO QCD with shower Monte Carlo algorithms”,  JHEP 11 (2004) 040, 10.1088/1126-6708/2004/11/040, arXiv:hep-ph/0409146.
  * [64] S. Frixione, P. Nason, and C. Oleari, “Matching NLO QCD computations with parton shower simulations: The powheg method”,  JHEP 11 (2007) 070, 10.1088/1126-6708/2007/11/070, arXiv:0709.2092.
  * [65] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: The POWHEG BOX”,  JHEP 06 (2010) 043, 10.1007/JHEP06(2010)043, arXiv:1002.2581.
  * [66] T. Melia, P. Nason, R. Rontsch, and G. Zanderighi, “$\PW^{+}\PW^{-}$, $\PW\PZ$ and $\PZ\PZ$ production in the POWHEG BOX”,  JHEP 11 (2011) 078, 10.1007/JHEP11(2011)078, arXiv:1107.5051.
  * [67] P. Nason and G. Zanderighi, “$\PW^{+}\PW^{-}$, $\PW\PZ$ and $\PZ\PZ$ production in the POWHEG-BOX-V2”,  Eur. Phys. J. C 74 (2014) 2702, 10.1140/epjc/s10052-013-2702-5, arXiv:1311.1365.
  * [68] J. M. Campbell, R. K. Ellis, P. Nason, and E. Re, “Top-pair production and decay at NLO matched with parton showers”,  JHEP 04 (2015) 114, 10.1007/JHEP04(2015)114, arXiv:1412.1828.
  * [69] S. Alioli, S.-O. Moch, and P. Uwer, “Hadronic top-quark pair-production with one jet and parton showering”,  JHEP 01 (2012) 137, 10.1007/JHEP01(2012)137, arXiv:1110.5251.
  * [70] S. Alioli, P. Nason, C. Oleari, and E. Re, “NLO single-top production matched with shower in POWHEG: $s$\- and $t$-channel contributions”,  JHEP 09 (2009) 111, 10.1088/1126-6708/2009/09/111, arXiv:0907.4076. [Erratum: 10.1007/JHEP02(2010)011].
  * [71] E. Re, “Single-top $Wt$ production matched with parton showers using the POWHEG method”,  Eur. Phys. J. C 71 (2011) 1547, 10.1140/epjc/s10052-011-1547-z, arXiv:1009.2450.
  * [72] R. Frederix, E. Re, and P. Torrielli, “Single-top $t$-channel hadroproduction in the four-flavour scheme with POWHEG and aMC@NLO”,  JHEP 09 (2012) 130, 10.1007/JHEP09(2012)130, arXiv:1207.5391.
  * [73] S. Alioli, P. Nason, C. Oleari, and E. Re, “NLO Higgs boson production via gluon fusion matched with shower in POWHEG”,  JHEP 04 (2009) 002, 10.1088/1126-6708/2009/04/002, arXiv:0812.0578.
  * [74] P. Nason and C. Oleari, “NLO Higgs boson production via vector-boson fusion matched with shower in POWHEG”,  JHEP 02 (2010) 037, 10.1007/JHEP02(2010)037, arXiv:0911.5299.
  * [75] K. Mimasu, V. Sanz, and C. Williams, “Higher order QCD predictions for associated Higgs production with anomalous couplings to gauge bosons”,  JHEP 08 (2016) 039, 10.1007/JHEP08(2016)039, arXiv:1512.02572.
  * [76] G. Heinrich et al., “NLO predictions for Higgs boson pair production with full top quark mass dependence matched to parton showers”,  JHEP 08 (2017) 088, 10.1007/JHEP08(2017)088, arXiv:1703.09252.
  * [77] G. Heinrich et al., “Probing the trilinear Higgs boson coupling in di-Higgs production at NLO QCD including parton shower effects”,  JHEP 06 (2019) 066, 10.1007/JHEP06(2019)066, arXiv:1903.08137.
  * [78] J. Davies et al., “Double Higgs boson production at NLO: Combining the exact numerical result and high-energy expansion”,  JHEP 11 (2019) 024, 10.1007/JHEP11(2019)024, arXiv:1907.06408.
  * [79] G. Heinrich, S. P. Jones, M. Kerner, and L. Scyboz, “A non-linear EFT description of $\Pg\Pg\to\PH\PH$ at NLO interfaced to POWHEG”,  JHEP 10 (2020) 021, 10.1007/JHEP10(2020)021, arXiv:2006.16877.
  * [80] T. Sjöstrand et al., “An introduction to PYTHIA $8.2$”,  Comput. Phys. Commun. 191 (2015) 159, 10.1016/j.cpc.2015.01.024, arXiv:1410.3012.
  * [81] B. Hespel, D. Lopez-Val, and E. Vryonidou, “Higgs pair production via gluon fusion in the two-Higgs-doublet model”,  JHEP 09 (2014) 124, 10.1007/JHEP09(2014)124, arXiv:1407.0281.
  * [82] A. Carvalho, “Gravity particles from warped extra dimensions, predictions for LHC”, 3, 2014. arXiv:1404.0102.
  * [83] NNPDF Collaboration, “Parton distributions for the LHC Run II”,  JHEP 04 (2015) 040, 10.1007/JHEP04(2015)040, arXiv:1410.8849.
  * [84] NNPDF Collaboration, “Parton distributions from high-precision collider data”,  Eur. Phys. J. C 77 (2017) 663, 10.1140/epjc/s10052-017-5199-5, arXiv:1706.00428.
  * [85] J. Butterworth et al., “PDF4LHC recommendations for LHC Run II”,  J. Phys. G 43 (2016) 023001, 10.1088/0954-3899/43/2/023001, arXiv:1510.03865.
  * [86] J. Rojo et al., “The PDF4LHC report on PDFs and LHC data: Results from Run I and preparation for Run II”,  J. Phys. G 42 (2015) 103103, 10.1088/0954-3899/42/10/103103, arXiv:1507.00556.
  * [87] A. Accardi et al., “A critical appraisal and evaluation of modern PDFs”,  Eur. Phys. J. C 76 (2016) 471, 10.1140/epjc/s10052-016-4285-4, arXiv:1603.08906.
  * [88] CMS Collaboration, “Event generator tunes obtained from underlying event and multiparton scattering measurements”,  Eur. Phys. J. C 76 (2016) 155, 10.1140/epjc/s10052-016-3988-x, arXiv:1512.00815.
  * [89] CMS Collaboration, “Extraction and validation of a new set of CMS PYTHIA$8$ tunes from underlying-event measurements”,  Eur. Phys. J. C. 80 (2020) 4, 10.1140/epjc/s10052-019-7499-4, arXiv:1903.12179.
  * [90] P. Skands, S. Carrazza, and J. Rojo, “Tuning PYTHIA $8.1$: The Monash 2013 tune”,  Eur. Phys. J. C 74 (2014) 3024, 10.1140/epjc/s10052-014-3024-y, arXiv:1404.5630.
  * [91] J. Alwall et al., “Comparative study of various algorithms for the merging of parton showers and matrix elements in hadronic collisions”,  Eur. Phys. J. C 53 (2008) 473, 10.1140/epjc/s10052-007-0490-5, arXiv:0706.2569.
  * [92] R. Frederix and S. Frixione, “Merging meets matching in MC@NLO”,  JHEP 12 (2012) 061, 10.1007/JHEP12(2012)061, arXiv:1209.6215.
  * [93] GEANT4 Collaboration, “: A simulation toolkit”,  Nucl. Instrum. Meth. A 506 (2003) 250, 10.1016/S0168-9002(03)01368-8.
  * [94] CMS Collaboration, “Particle-flow reconstruction and global event description with the CMS detector”,  JINST 12 (2017) P10003, 10.1088/1748-0221/12/10/P10003, arXiv:1706.04965.
  * [95] M. Cacciari, G. P. Salam, and G. Soyez, “The anti-$\kt$ jet clustering algorithm”,  JHEP 04 (2008) 063, 10.1088/1126-6708/2008/04/063, arXiv:0802.1189.
  * [96] M. Cacciari, G. P. Salam, and G. Soyez, “user manual”,  Eur. Phys. J. C 72 (2012) 1896, 10.1140/epjc/s10052-012-1896-2, arXiv:1111.6097.
  * [97] CMS Collaboration, “Electron and photon reconstruction and identification with the CMS experiment at the CERN LHC”,  JINST 16 (2021), no. 05, P05014, 10.1088/1748-0221/16/05/P05014, arXiv:2012.06888.
  * [98] CMS Collaboration, “Measurement of the Higgs boson production rate in association with top quarks in final states with electrons, muons, and hadronically decaying tau leptons at $\sqrt{s}=13\TeV$”,  Eur. Phys. J. C 81 (2021) 378, 10.1140/epjc/s10052-021-09014-x, arXiv:2011.03652.
  * [99] CMS Collaboration, “Performance of electron reconstruction and selection with the CMS detector in $\Pp\Pp$ collisions at $\sqrt{s}=8\TeV$”,  JINST 10 (2015) P06005, 10.1088/1748-0221/10/06/P06005, arXiv:1502.02701.
  * [100] CMS Collaboration, “Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at $\sqrt{s}=13\TeV$”,  JINST 13 (2018) P06015, 10.1088/1748-0221/13/06/P06015, arXiv:1804.04528.
  * [101] CMS Collaboration, “Performance of reconstruction and identification of $\PGt$ leptons decaying to hadrons and $\PGn_{\PGt}$ in $\Pp\Pp$ collisions at $\sqrt{s}=13\TeV$”,  JINST 13 (2018) P10005, 10.1088/1748-0221/13/10/P10005, arXiv:1809.02816.
  * [102] CMS Collaboration, “Identification of hadronic tau lepton decays using a deep neural network”, 1, 2022. arXiv:2201.08458. Submitted to JINST.
  * [103] Y. LeCun, “Generalization and network design strategies”, in  Connectionism in perspective, R. Pfeifer, Z. Schreter, F. Fogelman, and L. Steels, eds.  Elsevier, Zurich, Switzerland, 1989.  An extended version was published as a technical report of the University of Toronto.
  * [104] CMS Collaboration, “Pileup mitigation at CMS in $13\TeV$ data”,  JINST 15 (2020) P09018, 10.1088/1748-0221/15/09/p09018, arXiv:2003.00503.
  * [105] D. Bertolini, P. Harris, M. Low, and N. Tran, “Pileup per particle identification”,  JHEP 10 (2014) 059, 10.1007/JHEP10(2014)059, arXiv:1407.6013.
  * [106] CMS Collaboration, “Jet energy scale and resolution in the CMS experiment in $\Pp\Pp$ collisions at $8\TeV$”,  JINST 12 (2017) P02014, 10.1088/1748-0221/12/02/P02014, arXiv:1607.03663.
  * [107] CMS Collaboration, “Search for resonances decaying to a pair of Higgs bosons in the $\PQb\PAQb\PQq\PAQq^{\prime}\ell\PGn$ final state in proton-proton collisions at $\sqrt{s}=13\TeV$”,  JHEP 10 (2019) 125, 10.1007/JHEP10(2019)125, arXiv:1904.04193.
  * [108] E. Bols et al., “Jet flavour classification using DeepJet”,  JINST 15 (2020), no. 12, P12012, 10.1088/1748-0221/15/12/P12012, arXiv:2008.10519.
  * [109] CMS Collaboration, “Performance of missing transverse momentum reconstruction in proton-proton collisions at $\sqrt{s}=13\TeV$ using the CMS detector”,  JINST 14 (2019) P07004, 10.1088/1748-0221/14/07/P07004, arXiv:1903.06078.
  * [110] CMS Collaboration, “Evidence for associated production of a Higgs boson with a top quark pair in final states with electrons, muons, and hadronically decaying $\Pgt$ leptons at $\sqrt{s}=13\TeV$”,  JHEP 08 (2018) 066, 10.1007/JHEP08(2018)066, arXiv:1803.05485.
  * [111] L. Breiman, J. Friedman, R. A. Olshen, and C. J. Stone, “Classification and regression trees”.  Chapman and Hall/CRC, 1984.  ISBN 978-0-412-04841-8.
  * [112] K. Ehatäht, L. Marzola, and C. Veelken, “Reconstruction of the mass of Higgs boson pairs in events with Higgs boson pairs decaying into four $\PGt$ leptons”, 9, 2018. arXiv:1809.06140.
  * [113] T. Chen and C. Guestrin, “XGBoost: A scalable tree boosting system”, in  Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, p. 785.  ACM, New York, NY, USA, 2016.  arXiv:1603.02754.  10.1145/2939672.2939785.
  * [114] F. Pedregosa et al., “Scikit-learn: Machine learning in Python”, J. Mach. Learn. Res. 12 (2011) 2825, arXiv:1201.0490.
  * [115] L. Tani, D. Rand, C. Veelken, and M. Kadastik, “Evolutionary algorithms for hyperparameter optimization in machine learning for application in high energy physics”,  Eur. Phys. J. C 81 (2021) 170, 10.1140/epjc/s10052-021-08950-y, arXiv:2011.04434.
  * [116] LHC Higgs Cross Section Working Group, “Handbook of LHC Higgs Cross Sections: 4. Deciphering the nature of the Higgs sector”,  CERN Yellow Rep. Monogr. 2 (2017) 10.23731/CYRM-2017-002, arXiv:1610.07922.
  * [117] M. Grazzini, S. Kallweit, D. Rathlev, and M. Wiesemann, “$\PW^{\pm}\PZ$ production at hadron colliders in NNLO QCD”,  Phys. Lett. B 761 (2016) 179, 10.1016/j.physletb.2016.08.017, arXiv:1604.08576.
  * [118] F. Cascioli et al., “$\PZ\PZ$ production at hadron colliders in NNLO QCD”,  Phys. Lett. B 735 (2014) 311, 10.1016/j.physletb.2014.06.056, arXiv:1405.2219.
  * [119] Caola, F. and Melnikov, K. and Rontsch, R. and Tancredi, L., “QCD corrections to $\PZ\PZ$ production in gluon fusion at the LHC”,  Phys. Rev. D 92 (2015), no. 9, 094028, 10.1103/PhysRevD.92.094028, arXiv:1509.06734.
  * [120] CMS Collaboration, “Jet algorithms performance in $13\TeV$ data”, CMS Physics Analysis Summary CMS-PAS-JME-16-003, 2017.
  * [121] CMS Collaboration, “Measurement of the inelastic proton-proton cross section at $\sqrt{s}=13\TeV$”,  JHEP 07 (2018) 161, 10.1007/JHEP07(2018)161, arXiv:1802.02613.
  * [122] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, “Asymptotic formulae for likelihood-based tests of new physics”,  Eur. Phys. J. C 71 (2011) 1554, 10.1140/epjc/s10052-011-1554-0, arXiv:1007.1727. [Erratum: 10.1140/epjc/s10052-013-2501-z].
  * [123] ATLAS and CMS Collaborations, and LHC Higgs Combination Group, “Procedure for the LHC Higgs boson search combination in summer 2011”, technical report, 2011.
  * [124] J. S. Conway, “Incorporating nuisance parameters in likelihoods for multisource spectra”, in  PHYSTAT 2011, p. 115.  2011\.  arXiv:1103.0354.  10.5170/CERN-2011-006.115.
  * [125] T. Junk, “Confidence level computation for combining searches with small statistics”,  Nucl. Instrum. Meth. A 434 (1999) 435, 10.1016/S0168-9002(99)00498-2, arXiv:hep-ex/9902006.
  * [126] A. L. Read, “Presentation of search results: The $\textrm{CL}_{s}$ technique”,  J. Phys. G 28 (2002) 2693, 10.1088/0954-3899/28/10/313.
  * [127] HEPData record for this analysis, 2022.  10.17182/hepdata.130795.