Observed profile likelihood scans of $\kappa_\lambda$. The scans are performed by varying only the coupling modifier of interest, while all other relevant coupling modifiers are fixed to unity.

Expected profile likelihood scans of $\kappa_\lambda$. The scans are performed by varying only the coupling modifier of interest, while all other relevant coupling modifiers are fixed to unity.

Observed profile likelihood scans of $\kappa_{2V}$. The scans are performed by varying only the coupling modifier of interest, while all other relevant coupling modifiers are fixed to unity.

Expected profile likelihood scans of $\kappa_{2V}$. The scans are performed by varying only the coupling modifier of interest, while all other relevant coupling modifiers are fixed to unity.

Confidence level contours at 68% (solid line) and 95% (dashed line) in the $(\kappa_\lambda, \kappa_{2V})$ parameter space, when all other coupling modifiers are fixed to their SM predictions. The corresponding expected contours are shown by the inner and outer shaded regions The SM prediction is indicated by the star, while the best-fit value is denoted by the cross.

Confidence level contours at 68% (solid line) and 95% (dashed line) in the $(\kappa_\lambda, \kappa_{2V})$ parameter space, when all other coupling modifiers are fixed to their SM predictions. The corresponding expected contours are shown by the inner and outer shaded regions The SM prediction is indicated by the star, while the best-fit value is denoted by the cross.

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 2. 3he lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the high-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 2. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed $m_{\gamma\gamma}$ distributions for the low-mass categories in Run 3. The lines show the fit results for the continuum background only (light dotted), adding the single Higgs boson background (black dotted) and the full fit (solid).

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

Observed and expected profiled likelihood scans on $\kappa_\lambda$ for the low-mass (LM) and high-mass (HM) regions separately and for their combination. The scan is performed by floating the indicated parameter of interest while fixing all the others to their SM value.

The expected number of events (estimated by using simulation) from the SM HH signals and single Higgs boson production, and the expected number of events from the continuum background, evaluated in the 120 GeV < $m_{\gamma\gamma}$ < 130 GeV window in Run 3 for the different low-mass (LM) and high-mass (HM) categories. For comparison, the number of data events is also shown. The uncertainties in the HH signals and single Higgs boson backgrounds include the systematic uncertainties discussed in Section 6. Asymmetric uncertainties arise primarily from the theory calculation of the SM ggF HH cross section and the large uncertainty in the yield of single Higgs bosons produced in ggF events in association with heavy-flavour jets. The uncertainty in the continuum background is given by the sum in quadrature of the statistical uncertainty from the fit to the data and the spurious signal uncertainty.

The BDT distribution (with x-axis zoomed in) of the di-Higgs ggF signal for two different values of $\kappa_{\lambda}$ and the main backgrounds in the high mass region ($m^{*}_{b\bar{b}\gamma\gamma} > 350$ GeV). Distributions are normalized to unit area. The dotted lines denote the category boundaries. Events with a BDT score below 0.857 in the high mass region are discarded. The range of BDT scores is from 0.8 to 1.

The $BDT_{tot}$ score for the benchmark signal $m_{X}$ = 300 GeV and for the main backgrounds. Distributions are normalized to unit area. The dotted line denotes the event selection threshold. Events with a $BDT_{tot}$ score below 0.85 for $m_{X}$ = 300 GeV are discarded.

The $BDT_{tot}$ score for the benchmark signal $m_{X}$ = 500 GeV and for the main backgrounds. Distributions are normalized to unit area. The dotted line denotes the event selection threshold. Events with a $BDT_{tot}$ score below 0.75 for $m_{X}$ = 500 GeV are discarded.

Distributions of $m_{\gamma\gamma}$ in high mass BDT tight category for the nonresonant $HH$ search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband.

Distributions of $m_{\gamma\gamma}$ in high mass BDT loose category for the nonresonant $HH$ search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband.

Distributions of $m_{\gamma\gamma}$ in low mass BDT tight category for the nonresonant $HH$ search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband.

Distributions of $m_{\gamma\gamma}$ in low mass BDT loose category for the nonresonant $HH$ search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband.

Distributions of $m_{\gamma\gamma}$ for the selections used for the resonance mass point $m_{X}$ = 300 GeV for the resonant search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband. The scalar resonance signal is scaled to a total production cross section $\sigma(pp \rightarrow X \rightarrow HH)$ = 370 fb.

Distributions of $m_{\gamma\gamma}$ for the selections used for the resonance mass point $m_{X}$ = 500 GeV for the resonant search. The data-derived fractions of nonresonant $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and dijet background are applied and the total background is normalized to the data sideband. The scalar resonance signal is scaled to a total production cross section $\sigma(pp \rightarrow X \rightarrow HH)$ = 67 fb.

The number of data events observed in the 120 GeV < $m_{\gamma\gamma}$ < 130 GeV window, the number of $HH$ signal events expected for $\kappa_{\lambda}$ = 1 and for $\kappa_{\lambda}$ = 10, and events expected for single Higgs boson production (estimated using MC simulation), as well as for continuum background. For the single Higgs boson, Rest includes VBF, $WH$, $tHqb$, and $tHW$. The values are obtained from a fit of the Asimov data set generated under the SM signal-plus-background hypothesis, $\kappa_{\lambda}$ = 1. The continuum background component of the Asimov data set is obtained from the fit of the data sideband. The uncertainties in $HH$ signals and single Higgs boson background include the systematic uncertainties discussed in Section 6. The uncertainty in the continuum background is given by the sum in quadrature of the statistical uncertainty from the fit to the data and the spurious-signal uncertainty.

Observed and expected limits at 95% CL on the cross section of nonresonant Higgs boson pair production as a function of the Higgs boson self-coupling modifier $\kappa_{\lambda}= \lambda_{HHH}/\lambda^{\textrm{SM}}_{HHH}$. The expected constraints on $\kappa_{\lambda}$ are obtained with a background hypothesis excluding $pp \rightarrow HH$ production. The $\pm 1\sigma$ and $\pm 2\sigma$ variations about the expected limit due to statistical and systematic uncertainties are also shown. The theory prediction curve represents the scenario where all parameters and couplings are set to their SM values except for $\kappa_{\lambda}$. The uncertainty band of the theory prediction curve shows the cross-section uncertainty.

Values of the negative log-profile-likelihood ratio ($-2ln\Lambda$) as a function of $\kappa_{\lambda}$ evaluated for the combination of all the categories of the nonresonant search. The coupling of the Higgs boson to fermions and gauge bosons is set to SM values in the profile likelihood calculation. The Asimov data set is generated under the SM signal-plus-background hypothesis, $\kappa_{\lambda}$= 1. All systematic uncertainties, including the theoretical uncertainties in the di-Higgs boson production cross section, are included. The intersections of the solid curves and the horizontal dashed lines indicate the 1$\sigma$ and 2$\sigma$ confidence-level intervals.

The number of events observed in the 120 < $m_{\gamma\gamma}$ < 130 GeV window in data, the number of events expected for scalar resonance signals of masses $m_{X}$ = 300 GeV and $m_{X}$ = 500 GeV assuming a total production cross section $\sigma(pp \rightarrow X \rightarrow HH)$ equal to the observed exclusion limits of Figure 15, and events expected for SM $HH$ and single Higgs boson production (estimated using MC simulation), as well as for continuum background. The values are obtained from a fit of the Asimov data set generated under the signal-plus-background hypothesis. The continuum background component of the Asimov data set is obtained from the fit of the data sideband. The uncertainties in the resonant signals and the SM $HH$ and single-Higgs-boson backgrounds include the systematic uncertainties discussed in Section 6. The uncertainty in the continuum background is given by the sum in quadrature of the statistical uncertainty from the fit to the data and the spurious-signal uncertainty.

Observed and expected limits at 95% CL on the production cross section of a narrow-width scalar resonance $X$ as a function of the mass $m_{X}$ of the hypothetical scalar particle. The black solid line represents the observed upper limits. The dashed line represents the expected upper limits. The $\pm 1\sigma$ and $\pm 2\sigma$ variations about the expected limit due to statistical and systematic uncertainties are also shown.

Breakdown of the dominant systematic uncertainties. The impact of the uncertainties is defined according to the statistical analysis described in Section 7. It corresponds to the relative variation of the expected upper limit on the cross section when re-evaluating the profile likelihood ratio after fixing the nuisance parameter in question to its best-fit value, while all remaining nuisance parameters remain free to float. The impact is shown in %. Only systematic uncertainties with an impact of at least 0.2% are shown. Uncertainties of the "Norm. + Shape" type affect both the normalization and the parameters of the functional form. The rest of the uncertainties affect only the yields.

Cutflow for nonresonant di-Higgs ggF signal sample, yields are normalized to 139 $fb^{-1}$.

Cutflow for resonant signal sample, with $m_{X}$ = 300 GeV, yields are normalized to 139 $fb^{-1}$.

Cutflow for resonant signal sample, with $m_{X}$ = 500 GeV, yields are normalized to 139 $fb^{-1}$.

Comparison of $m_{b\bar{b}}$ distributions when applying the specific b-jet energy calibration and the nominal jet energy calibration. The distributions are fitted using a Bukin function, and the values of the peak position, resolution and the relative improvement are reported in the legend.

Fit results of $m_{b\bar{b}}$ distributions when applying the specific b-jet energy calibration and the nominal jet energy calibration. The distributions are fitted using a Bukin function, and the values of the peak position, resolution and the relative improvement are reported in the legend.

The relative amount (purity) of expected events from SM $HH$ and single Higgs boson production processes for each of the four categories of the nonresonant search. The Higgs boson pair production with $k_{\lambda} = 1$ is considered as signal in (a), while the case with $k_{\lambda} = 10$ is considered as signal in (b).

The expected significance in each of the four categories of the nonresonant search. The Higgs boson pair production with $k_{\lambda} = 1$ is considered as signal in (a), while the case with $k_{\lambda} = 10$ is considered as signal in (b). The single Higgs boson processes and the di-photon continuum spectrum are considered as background.

Spurious signal result for the exponential pdf for the various ggF nonresonant di-Higgs categories. In each category, the spurious signal value ($n_{sp}$) and its ratio to the expected statistical error ($Z_{spur}$) from data are shown.

Spurious signal result for the exponential pdf as function of the resonant di-Higgs signal mass. The spurious signal value and its ratio to the expected statistical error from data are shown.

Distributions of the signal efficiency as a function of $\kappa_{\lambda}$, for the di-Higgs boson ggF nonresonant production mode. The range of $\kappa_{\lambda}$ in the table is from -10 to 10.

Distributions of the signal efficiency as a function of $\kappa_{\lambda}$, for the di-Higgs boson VBF nonresonant production mode. The range of $\kappa_{\lambda}$ in the table is from -10 to 10.

Values of the negative log-profile-likelihood ($-2ln\Lambda$) as a function of $\kappa_{\lambda}$ evaluated for the combination of all the categories of the nonresonant search. The coupling of the Higgs boson to fermions and gauge bosons is set to SM values in the profile likelihood calculation. The Asimov data set is generated under the SM signal-plus-background hypothesis, $\kappa_{\lambda}$= 1. All systematic uncertainties, including the theoretical uncertainties on the di-Higgs boson production cross section, are included. The intersections of the solid curves and the horizontal dashed lines indicate the 1$\sigma$ and 2$\sigma$ confidence level intervals. The best fit value corresponds to $\kappa_{\lambda}$ = $2.8^{+2.0}_{-2.2}(^{+3.8}_{-4.3})$ for the 1$\sigma$(2$\sigma$) confidence interval. The expected value corresponds to $\kappa_{\lambda}$ = $1.0^{+5.5}_{-2.4}(^{+7.3}_{-4.2})$ for the 1$\sigma$(2$\sigma$) confidence interval. The dashed curves represent values of the negative log-profile-likelihood where the Higgs boson branching fractions and the cross section of the production modes are varied as a function of $\kappa_{\lambda}$. In this case,the best fit value corresponds to $\kappa_{\lambda}$ = $2.7^{+2.0}_{-2.2}(^{+3.8}_{-4.3})$ and the expected value corresponds to $\kappa_{\lambda}$ = $1.0^{+5.4}_{-2.5}(^{+7.3}_{-4.3})$ for the 1$\sigma$(2$\sigma$) confidence interval.

Minimum BDT value of the events passing the selection criteria of the resonant search. The combined BDT score is formed using as coefficients $C_{1}$ = 0.65 and $C_{2}$ = 1 − $C_{1}$. The selection efficiency for the resonant $X \rightarrow HH$ signal is also shown.