Double Higgs boson production and Higgs self-coupling extraction at CLIC

Abstract

The Compact Linear Collider (CLIC) is a future electron–positron collider that will allow measurements of the trilinear Higgs self-coupling in double Higgs boson events produced at its high-energy stages with collision energies from $\sqrt{s}$ = 1.4 to 3 TeV. The sensitivity to the Higgs self-coupling is driven by the measurements of the cross section and the invariant mass distribution of the Higgs-boson pair in the W-boson fusion process, $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$. It is enhanced by including the cross-section measurement of ZHH production at 1.4 TeV. The expected sensitivity of CLIC for Higgs pair production through W-boson fusion is studied for the decay channels $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  using full detector simulation including all relevant backgrounds at $\sqrt{s}$ = 1.4 TeV with an integrated luminosity of $\mathcal {L}$ = 2.5 ab$^{-1}$ and at $\sqrt{s}$ = 3 TeV with $\mathcal {L}$ = 5 ab$^{-1}$. Combining $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ and ZHH  cross-section measurements at 1.4 TeV with differential measurements in $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ events at 3 TeV, CLIC will be able to measure the trilinear Higgs self-coupling with a relative uncertainty of $-8\%$ and $+11\%$ at 68% C.L., assuming the Standard Model. In addition, prospects for simultaneous constraints on the trilinear Higgs self-coupling and the Higgs-gauge coupling HHWW are derived based on the ${\text {H}\text {H}\nu \bar{\nu }}$ measurement.

1 Introduction

The discovery of the Higgs boson [1, 2] has initiated an era of investigations of its properties and of the nature of the mechanism that breaks the electroweak symmetry. Besides its mass and width, the properties of interest include the couplings of the Higgs boson to other Standard Model (SM) and hypothetical non-SM particles as well as the coupling to itself. While the couplings to other SM particles illustrate the way these particles obtain masses in the Higgs mechanism, the self-coupling parameter determines the shape of the Higgs potential which has implications for the vacuum metastability, the hierarchy problem, as well as the electroweak phase transition and baryogenesis. In the Standard Model, the Higgs potential for the Higgs field $\phi$ is described by

$$\begin{aligned} V(\phi ) = \mu ^{2} \phi ^{\dagger } \phi + \frac{\lambda ^2}{2} (\phi ^{\dagger } \phi )^{2}, \end{aligned}$$

(1)

where $\mu$ is proportional to the Higgs boson mass and $\lambda$ is the Higgs self-coupling. This implies a fixed relation $m_H^2 = \lambda v$ between the mass and the self-coupling, with the vacuum expectation value v. In the interaction Lagrangian, this potential leads to a trilinear self-coupling ${g_{\text {H}\text {H}\text {H}}}$ which is proportional to $\lambda$.

A deviation of the Higgs potential from the SM would directly point to new physics, for example in the context of baryogenesis: Indeed, one of the conditions for electroweak baryogenesis is the presence of a strong first-order phase transition in the breaking of the electroweak symmetry in the early universe. In order to modify the Higgs potential accordingly, at least one additional scalar needs to be introduced [3]. This can be an additional scalar singlet [4] or doublet. The latter is realised in two-Higgs doublet models (2HDM) [5, 6], which introduce four additional scalars. These models can lead to modifications of the Higgs self-coupling. A sufficiently heavy neutral scalar can cause a resonance in the invariant mass of the Higgs boson pair in the production of two SM-like Higgs bosons.

Existing models including those discussed above with additional scalars as well as theories where the Higgs boson is composite predict differences of the Higgs self-couplings to the SM value between a few and tens of percent [7]. These estimates assume the scenario that no additional states of the electroweak symmetry breaking sector can be discovered at the LHC. An overview of BSM theories modifying the Higgs self-coupling is given in [8].

A measurement of the Higgs self-coupling with a precision of better than 50% will not be possible at the High-Luminosity Large Hadron Collider (HL-LHC) [9]. More precise measurements are possible at high-energy linear colliders, as they give direct access to double Higgs boson production in a comparably clean environment. Electron–positron colliders below a center-of-mass energy of $\sqrt{s} \approx 500$ GeV do not have access to double Higgs boson production. They can only constrain the Higgs self-coupling indirectly through its loop contributions to single Higgs boson production [10]. The prospects for several proposed future options are discussed in [11]. The potential of the International Linear Collider (ILC) to measure the Higgs self-coupling directly in double Higgs boson production in association with a Z boson at $\sqrt{s}=500$ GeV and in the W-boson fusion double Higgs production channel at 1 TeV is described in [12,13,14].

The Compact Linear Collider (CLIC) is a mature option for a future linear electron–positron collider [15], which will allow the precise determination of the properties of the Higgs boson well beyond the precision of the HL-LHC. A detailed investigation of the CLIC prospects for the Higgs couplings to SM particles is given in [16] and an update of these results to a new luminosity and polarisation baseline scenario is provided in [17, 18]. A preliminary study of the Higgs self-coupling measurement at CLIC, based only on the measurement of the double Higgs boson production, has been presented in [16]. The analysis is updated and extended in this paper, most importantly by exploiting differential distributions in the analysis of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ at 3 TeV, by illustrating the impact of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {Z}\text {H}\text {H}}$ at 1.4 TeV, and by extracting ${g_{\text {H}\text {H}\text {H}}}$ and ${g_{\text {H}\text {H}\text {W}\text {W}}}$ in a joint fit.

Each energy stage at CLIC contributes to the indirect measurement of the Higgs self-coupling in single Higgs boson production. Combined with the HL-LHC standalone precision of 47% in the one-parameter fit, the CLIC run at the collision energy of 380 GeV will only improve this precision to 46% [11]. With increasing statistics and energy, the indirect limits in the one-parameter fit will be improved to 41% after the energy stage at $\sqrt{s}$ = 1.5 TeV and 35% after the 3 TeV energy stage. However, already at 1.5 TeV, the direct accessibility of double Higgs production allows much more powerful, potentially model-independent constraints to be put on the Higgs self-coupling which by far exceed the precision obtained in single Higgs measurements [11]. These measurements are the subject of this paper.

The high-energy stages of CLIC with centre-of-mass energies of 1.5 and 3 TeV provide the opportunity to access directly the trilinear Higgs self-coupling in double Higgs boson production. In the present study, the earlier choice of centre-of-mass energy for the second stage of 1.4 TeV [19] is used due to the availability of full simulation event samples. While we therefore base the following study on a run at 1.4 TeV, the prospects for 1.5 TeV are expected be very similar. The main channels are double Higgsstrahlung ZHH  production at 1.4 TeV and double Higgs boson production via W-boson fusion at 1.4 and 3 TeV. Both are directly sensitive to the trilinear Higgs self-coupling ${g_{\text {H}\text {H}\text {H}}}$, while the latter is also sensitive to the quartic Higgs-gauge coupling ${g_{\text {H}\text {H}\text {W}\text {W}}}$. This paper uses full detector simulation to study the CLIC potential for extracting these couplings from measurements of double Higgs boson production.

In a full Effective Field Theory approach, other operators apart from the one modifying the triple Higgs vertex can also contribute to the same final state. As these operators are themselves constrained by other measurements, e.g. single Higgs boson production channels, a global fit approach as studied in [20] is appropriate. Results for CLIC are presented in [21, Sec. 2.2.1], showing that the constraints from the global fit are very close to the ones obtained in the exclusive approach, due to the high precision measurements of other processes at CLIC. A detailed study of the impact of other operators was performed for the ZHH channel at 500 GeV [22].

This paper investigates the prospects for extracting the trilinear Higgs self-coupling at CLIC in double Higgs boson production at the high-energy stages of CLIC. It is structured as follows: Sect. 2 describes the strategy of the analysis and the various contributions to the sensitivity. In Sect. 3, the definition of the signal and background processes, as well as the simulation and reconstruction chain, are described. The event selection procedures for the analyses at 1.4 and 3 TeV for ${\text {H}\text {H}\nu \bar{\nu }}$ $\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ $\nu \bar{\nu }$ and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ $\nu \bar{\nu }$ are explained in Sect. 4. This is followed by the results for the cross section measurement in Sect. 5 and for the differential measurement giving the most stringent constraints in Sect. 6. A summary is provided in Sect. 7.

2 Analysis strategy

At CLIC, the Higgs self-coupling can be directly accessed through the measurement of double Higgs boson production. Two main channels contribute: W-boson fusion (WBF) double Higgs boson production (dominant part of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$) and the double Higgsstrahlung process ($\text {e}^{+}\text {e}^{-}\rightarrow {\text {Z}\text {H}\text {H}}$). The other process of vector boson fusion, namely Z-boson fusion ($\text {e}^{+}\text {e}^{-}\rightarrow \text {H}\text {H}\text {e}^{+}\text {e}^{-}$), has a one order of magnitude smaller cross section and is therefore not considered here. The dependence of the cross section on the centre-of-mass energy obtained with Whizard 1.95 [23, 24] is shown in Fig. 1. This illustrates that the highest cross section of ZHH  production among the forseen CLIC energy stages is at the first one above 500 GeV, assumed to be at 1.4 TeV in this paper. In ${\text {Z}\text {H}\text {H}}$  production, this energy stage also gives the best sensitivity to ${g_{\text {H}\text {H}\text {H}}}$. The cross section of WBF double Higgs boson production grows with the collision energy. Therefore, assuming the same polarisation configuration, the 3 TeV stage gives the largest event rate of WBF double Higgs boson production at CLIC. In $\text {e}^{+}\text {e}^{-}$ collisions at $\sqrt{s} \gtrsim 1.2$ TeV, WBF is the dominant double Higgs boson production mode for unpolarised beams. Its total cross section at 3 TeV, including effects of the luminosity spectrum and initial state radiation, exceeds that of double Higgsstrahlung at 1.4 TeV by a factor of 6. The single most sensitive measurement of Higgs boson pair production at CLIC is therefore the double Higgs boson production through WBF at 3 TeV.

Fig. 1

Cross section as a function of centre-of-mass energy for $\text {e}^{+}\text {e}^{-}\rightarrow {\text {Z}\text {H}\text {H}}$ and $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ production for a Higgs boson mass of $m_{\mathrm {H}}$ = 126 GeV. The values shown correspond to unpolarised beams including initial state radiation but not including the effect of beamstrahlung [16]

Fig. 2

Main Feynman diagrams contributing to double Higgs boson production via W-boson fusion. Diagram a contains the trilinear Higgs self-coupling, b grows with the quartic coupling ${g_{\text {H}\text {H}\text {W}\text {W}}}$, while c, d are sensitive to the Higgs coupling to W bosons

Figure 2 shows the main Feynman diagrams contributing to Higgs pair production via W-boson fusion. This channel contains the ${\text {H}\text {H}\text {H}}$ vertex which depends on the trilinear Higgs self-coupling ${g_{\text {H}\text {H}\text {H}}}$, as well as the ${\text {H}\text {H}\text {W}\text {W}}$ vertex which depends on the quartic Higgs-gauge coupling ${g_{\text {H}\text {H}\text {W}\text {W}}}$. Deviations from the SM values are defined as:

$$\begin{aligned} {\kappa _{\text {H}\text {H}\text {H}}}:= \frac{{g_{\text {H}\text {H}\text {H}}}}{{g_{\text {H}\text {H}\text {H}}^{\text {SM}}}} \mathrm {~and~} {\kappa _{\text {H}\text {H}\text {W}\text {W}}}:= \frac{{g_{\text {H}\text {H}\text {W}\text {W}}}}{{g_{\text {H}\text {H}\text {W}\text {W}}^{\text {SM}}}}. \end{aligned}$$

The total cross sections of WBF and Higgsstrahlung double Higgs boson production are sensitive to the value of the trilinear Higgs self-coupling. Figure 3 shows the parabolic dependence of the WBF double Higgs boson production cross section on ${g_{\text {H}\text {H}\text {H}}}$ at 3 TeV. The cross section at around $2.3\times {g_{\text {H}\text {H}\text {H}}^{\text {SM}}}$ is identical to the SM cross section. Therefore, only measuring the total cross section of this process will not be sufficient to determine ${g_{\text {H}\text {H}\text {H}}}$ unambiguously. This can be resolved by measuring the double Higgsstrahlung cross section which has an unambiguous dependence on ${g_{\text {H}\text {H}\text {H}}}$ as illustrated in Fig. 3 for 1.4 TeV. Another way to resolve the ambiguity is by using differential distributions such as the di-Higgs invariant mass [25]. It can also be exploited to distinguish whether a possible deviation from the SM originates from a modification of the ${\text {H}\text {H}\text {H}}$ or of the ${\text {H}\text {H}\text {W}\text {W}}$ vertex [25]. Differential distributions are therefore used in the following analysis (Sect. 6).

This analysis is focused on the two decay channels $\text {H}\text {H}\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ (branching fraction 34%) and $\text {H}\text {H}\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ $\rightarrow \mathrm{b}\bar{\mathrm{b}}\mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$ (branching fraction 8.4%). Both channels benefit from the relatively clean environment in electron–positron collisions at CLIC, the excellent jet energy resolution of the assumed CLIC detector concept using particle flow analysis, as well as from its very good flavour tagging capabilities [15]. This allows reconstruction of the kinematic properties of the Higgs boson pair.

The baseline scenario for CLIC sets the collision energy of the second stage to 1.5 TeV [17]. The earlier choice of 1.4 TeV [19] is used in the present study. It is expected that prospects for 1.5 TeV will be very similar as the cross section only changes by $-7\%$ for ${\text {Z}\text {H}\text {H}}$  and $+18\%$ for ${\text {H}\text {H}\nu \bar{\nu }}$. Results presented here are based on an integrated luminosity of 2.5 ab$^{-1}$ at a centre-of-mass energy of 1.4 TeV and 5 ab$^{-1}$ at $\sqrt{s}$ = 3 TeV.

The CLIC electron beam can be polarised with a polarisation of up to ± 80%. The negative polarisation of $-80\%$ leads to an increase of the cross section for $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\mathrm {\nu }_{\text {e}}\bar{\nu }_{\text {e}}}$ by a factor of 1.8. The positive polarisation has the inverse effect of reducing the cross section to 20% [16]. For the process $\text {e}^{+}\text {e}^{-}\rightarrow$ ${\text {Z}\text {H}\text {H}}$, the cross-section scaling factors are 1.12 (0.88) for the electron beam polarisation of $-80\%$ (+ 80%). Running a fraction of the integrated luminosity with positive polarisation is, however, desirable for other measurements including two-fermion production [21]. Therefore, a scheme of collecting 80% (20%) of the data with $-80\%$ (+ 80%) electron beam polarisation is envisaged, which is denoted by “4:1 polarisation scheme” in the following. A polarisation scaling factor $f_p$ is defined as the ratio of the total number of events for the assumed polarisation running scheme with respect to the total number of events without beam polarisation for the same total luminosity. We apply these scaling factors to obtain the total number of signal and background events for the entire energy stage. The treatment of the polarisation is detailed in Sect. 5.1. A proper optimisation of the selection criteria taking into account the polarisation dependent kinematics would result in a better signal selection and hence a higher significance. The chosen approach is conservative compared to a proper combination of data sets.

Fig. 3

Cross section dependence of the trilinear Higgs self-coupling for the processes ${\text {H}\text {H}\nu \bar{\nu }}$ production at 1.4 and 3 TeV and ${\text {Z}\text {H}\text {H}}$  production at 1.4 TeV for unpolarised beams. Beamstrahlung and initial state radiation are included. The SM case is at $\frac{{g_{\text {H}\text {H}\text {H}}}}{{g_{\text {H}\text {H}\text {H}}^{\text {SM}}}} =1$. The ambiguity of the cross section value in the case of ${\text {H}\text {H}\nu \bar{\nu }}$ production is illustrated for 3 TeV. No such ambiguity exists in the ${\text {Z}\text {H}\text {H}}$  production process

3 Simulation and reconstruction for signal and background samples

3.1 Definition of signal and background processes and Monte Carlo generation

The process $\text {e}^{+}\text {e}^{-}$ $\rightarrow$ ${\text {H}\text {H}\nu \bar{\nu }}$ with a total cross section of 0.59 fb (0.149 fb) at $\sqrt{s}=$ 3 TeV (1.4 TeV) in the decay channels $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ defines the signal. This includes a contribution from $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$  which cannot be distinguished from WBF experimentally. It amounts to a fraction of 1.76% of the total $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ cross section for the unpolarised case at 3 TeV. In the baseline polarisation scheme at 1.4 TeV, this contribution is larger: 13.5% for unpolarised beams, 9.3% and 39% for negatively and positively polarised electron beams, respectively. However, these ratios are still small compared to the statistical uncertainty of the measurement, keeping in mind that the case of positively polarised beams, which has the largest contribution of $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$ to the ${\text {H}\text {H}\nu \bar{\nu }}$ final state, only contributes 20% of the luminosity collected at 1.4 TeV. The background consists of processes with multiple intermediate electroweak gauge bosons resulting in multiple jets, single Higgs boson production in association with electroweak gauge bosons decaying to hadrons, as well as di-Higgs production with decays to other final states. In order to avoid overlap, Higgs boson pair production is removed from the inclusive multi-quark background samples. Specifically, the background processes which turned out to be non-negligible after the selection are $\text {e}^{+}\text {e}^{-}\rightarrow \mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$ (only relevant at 3 TeV), $\text {e}^{+}\text {e}^{-}\rightarrow \mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}\nu \bar{\nu }$, $\text {e}^{+}\text {e}^{-}\rightarrow \mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}\text {l}\bar{\mathrm {\nu }}$, $\text {e}^{+}\text {e}^{-}\rightarrow \mathrm{q}\bar{\mathrm{q}}\text {H}\nu \bar{\nu }$, $\text {e}^{\pm } \upgamma \rightarrow \mathrm {\nu }\mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$, and $\text {e}^{\pm } \upgamma \rightarrow \mathrm{q}\bar{\mathrm{q}}\text {H}\mathrm {\nu }$ where $\text {q}$ refers to $\text {u}$, $\text {d}$, $\text {s}$, $\text {c}$, and $\text {b}$ quarks, $\text {l}$ = $\text {e}^\pm$, $\mathrm {\mu }^\pm$, $\mathrm {\tau }^\pm$, and $\mathrm {\nu }= \mathrm {\nu }_{\text {e}}, \mathrm {\nu }_{\mathrm {\mu }}, \mathrm {\nu }_{\mathrm {\tau }},$ as well as the respective anti-particles ($\bar{\mathrm{q}}$  and $\bar{\mathrm {\nu }}$). The processes which do not contain explicitly a Higgs boson in the final state do not include Higgs propagators. All SM Higgs boson decays were included otherwise.

Initial state radiation [26] and beamstrahlung [27, 28] lead to a tail in the distribution of the effective centre-of-mass energy, which is included in the simulation. In addition to $\text {e}^{+}\text {e}^{-}$ collisions, photon-initiated processes are also considered. Processes with photons from beamstrahlung in the initial state are normalised to the corresponding lower luminosity. “Quasi-real” photons are modeled using the Equivalent Photon Approximation [29,30,31] as implemented in Whizard 1.95 [23, 24].

The contributions of the most important background processes are presented in Tables 1, 2,  3 and 4.

All samples are generated with Whizard 1.95 interfaced to Pythia 6.4 [32] for parton shower and hadronisation as well as Higgs decays. Tauola [33, 34] is used for $\tau$ lepton decays.

Unpolarised beams are assumed in the simulation samples. We have studied the effects of polarisation on the kinematics of a process where a potentially relevant effect is expected: The $\text {e}^{+}\text {e}^{-}\rightarrow \text {W}\text {W}$ background strongly decreases with positive electron beam polarisation as the contribution of the t-channel neutrino exchange is suppressed. However, while the kinematic distributions differ between 100% positive (only s-channel diagrams) and 100% negative (s- and t-channel diagrams) beam polarisation, the contribution from negatively polarised electrons dominates by far in both the $P(e^{-})=-80\%$ and $P(e^{-})=+80\%$ beam polarisation modes. Therefore, the W boson kinematics are unchanged, such that only the different normalisation between positive and negative beam polarisation modes has been taken into account in this study.

3.2 Detector simulation

The simulation in this analysis uses the CLIC_ILD detector model [15]. It is based on the ILD detector concept [35, 36] for the International Linear Collider (ILC) [37] adapted to the experimental conditions at CLIC: Due to the higher collision energy at CLIC than at the ILC, jets tend to be more energetic. Therefore, the hadronic calorimeter has more interaction lengths in CLIC_ILD than in the ILD concept (7.5 instead of 5.5 $\uplambda _{\text {I}}$). The magnetic field is slightly higher (4 instead of 3.5 T). The inner radius of the vertex detector is 31 mm in CLIC_ILD and 16 mm in ILD. In addition, the very forward detectors for CLIC_ILD have been redesigned as the detector at CLIC is required to cope with more beam-induced background in particular in the forward region. The CLIC_ILD detector has a cylindrical layout. The innermost subdetector is an ultra-light silicon vertex detector with six layers with a single point resolution of 3 $\upmu$m. It is surrounded by a large tracking system consisting of a large central gaseous Time Projection Chamber (TPC) surrounded by several silicon strip layers. Highly granular electromagnetic and hadronic calorimeters are located around the tracker. They are optimised for particle flow analysis which aims at reconstructing the final-state particles within a jet using the information from the tracking detectors combined with that from the calorimeters. The outermost part of the detector consists of an iron return yoke, which is instrumented with muon chambers. The forward region is equipped with a system of two electromagnetic calorimeters, the BeamCal and LumiCal. They are specifically designed for the luminosity measurement and the identification of electromagnetic clusters from forward electrons or photons.

At CLIC, the bunch crossings are separated by 0.5 ns. At the 3 TeV stage, there are on average 3.2 $\mathrm {\gamma \gamma }$  $\rightarrow$ hadrons interactions per bunch crossing [15]. In order to suppress the beam-induced background collected over the duration of a bunch train, the hit time resolution in the calorimeters is 1 ns while the TPC integrates over the entire bunch train. The elements of the silicon envelope of the TPC and the vertex detector have a time resolution of $10/\sqrt{12}$ ns.

Recently, a new detector model, CLICdet, has been optimised and validated for CLIC [38]. The performance of this analysis is expected to be similar if the CLICdet model had been used.

The detector simulation of the generated event samples is performed with Geant4 [39, 40] and the detector description toolkit Mokka [41]. Hits from 60 bunch crossings of beam-induced $\mathrm {\gamma \gamma }$  $\rightarrow$ hadrons background are overlaid to each event. This is done for all subdetectors. For most of them, this is more than the reconstruction window and hit resolution requires. For the TPC, this is a compromise between realism and computing capacities [15, 42].

3.3 Reconstruction

The reconstruction algorithms run in the Marlin framework [43] which is a part of iLCSoft [44]. This includes track reconstruction with the ILD track reconstruction software [45] and particle flow analysis based on tracks and calorimeter deposits with the PandoraPFA program [46,47,48] resulting in Particle Flow Objects (PFOs). Cuts on the timing of the PFOs are applied to suppress beam-induced backgrounds from other bunch crossings. Muon and electron candidates are identified using calorimeter and tracking information. They are required to be isolated by applying quality criteria on their impact parameters and by restricting the energy in the surrounding cone in dependence on the track energy. As the forward calorimeters were not used in the reconstruction, the geometrical acceptance and the efficiency of the forward calorimeters BeamCal and LumiCal from dedicated full simulation [49] are used to simulate the veto of forward electrons occurring in background processes in the polar angle region between 10 and 110 mrad.

Jets are reconstructed using the FastJet [50] package via the MarlinFastJet interface. Both the VLC algorithm [51, 52]¹ and the longitudinally invariant $k_{t}$ algorithm [53] are used in the analysis. The parameter settings for the jet reconstruction in the individual channels are specified in Sect. 4. Vertex reconstruction and heavy-flavour tagging is performed using the Linear Collider Flavour Identification (LcfiPlus) program [54]. Hadronic tau decays are identified using the TauFinder package [55].

The jets studied in this paper are predominantly b-jets with an energy around 100 GeV which are rather forward in the detector. The pure relative jet energy resolution achievable with the CLIC_ILD detector is between 3 and 5% for light-flavour jets [15, 47]. However, for forward b-jets such as those in this analysis, the resolution is degraded for several reasons: A part of the jet energy is missing due to neutrinos from heavy flavour decays and due to forward particles outside of the detector acceptance. In addition, the beam-induced background is higher in the forward region.

The momentum resolution for forward tracks with a transverse momentum around 100 GeV is estimated to be around $\sigma (\Delta p_{\text {T}}/p_{\text {T}}^2) = 9\times 10^{-4}$ GeV$^{-1}$ [15] and the impact parameter resolution is $\sigma _{d_0} \approx 1.5\,\upmu$m for central tracks and around $\sigma _{d_0} \approx 3\,\upmu$m in the forward region [15]. The b-tagging performance is expected to provide a mis-identification rate of around 0.1% for light flavor jets with a b-tagging efficiency of 55% [56, Fig. 6] for jets with a similar polar angle distribution as the signal.

4 Event selection

4.1 Common preselection and definition of orthogonal samples

To select events originating from double Higgs production in the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ $\rightarrow \mathrm{b}\bar{\mathrm{b}}\mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$ decay channels, all events containing isolated leptons (electrons, muons or hadronic $\tau$ leptons) are rejected. For this, electron and muon candidates compatible with prompt production with an absolute impact parameter below 0.04 mm (0.06 mm) for electrons (muons) are used. Furthermore, the fraction of energy deposited in the electromagnetic calorimeter $R_{\text {cal}}$ is required to be $R_{\text {cal}} > 0.9$ for electrons and $0.05< R_{\text {cal}} < 0.25$ for muons. A minimum track energy of 15 GeV is required, and an energy-dependent cone-based isolation criterion is imposed, allowing a typical maximum energy in the cone of, e.g., 23 GeV for 100 GeV tracks. For the identification of hadronically decaying $\tau$ leptons, parameters are chosen to optimise the performance for this analysis. In particular, a maximum energy of 3 GeV in the cone between 0.03 and 0.33 rad around the seed particle is required.

In order to define orthogonal samples to be used for the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ channels, the events are clustered into four jets using the $k_{t}$ algorithm with a jet size parameter of $R=0.7$. A flavour tagging algorithm is applied on these jets using the LcfiPlus package. It first identifies the primary vertices, followed by the secondary vertices indicating b and c hadron decays. Then, the secondary vertices are assigned to jets. In the next step, the jet clustering is refined by using as seeds only those tracks and leptons originating from secondary vertices. Finally, values for b tags, c tags, and light-flavour quark tags are assigned to each jet. This classification is based on a multivariate discriminant trained on $\text {e}^{+}\text {e}^{-}\rightarrow \text {Z}\nu \bar{\nu }$ events, which have a similar event topology to the signal events. In the training, events with Z bosons subsequently decaying to $\mathrm{b}\bar{\mathrm{b}}$ are treated as signal, while those decaying to either $\mathrm{c}\bar{\mathrm{c}}$ or $\mathrm{q}\bar{\mathrm{q}}$ (with q = $\text {u}$, $\text {d}$, $\text {s}$) are considered background. The b-tagging performance relies on the ability to identify secondary vertices and tracks which do not originate from the primary interaction point. This depends in particular on the single point resolution of the vertex detector. In the CLIC_ILD model, this is assumed to be $\approx$ 3 $\upmu$m. Flavor tagging performances reached with the CLIC_ILD detector at 1.4 TeV are illustrated in [56, Fig. 6]. The $\sum _{4}{b\text {-tag}}$ distribution at 3 TeV is shown in Fig. 4. The $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  final state tends to be in the region between 2 and 4 of the $\Sigma _4$ $b$-tag distribution, which is much higher than the backgrounds. Contributions from other Higgs decays tend to values between 0 and 2.5. This shows that this criterion can be used to remove background contributions, and a large contribution of other $\text {H}\text {H}$ decay channels is also removed. The sample is then split into mutually exclusive samples with $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ candidates in the following way: Events are chosen as $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ candidates if the sum of the b-tag values $\sum _{4}{b\text {-tag}}$ of the jets is smaller than 1.5 (2.3) at 1.4 TeV (3 TeV). Otherwise, the events are considered as $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ candidates. Further selection criteria are applied separately for the two channels.

Fig. 4

Distribution of the sum of the b-tag values for the inclusive ${\text {H}\text {H}\nu \bar{\nu }}$ and the ${\text {H}\text {H}\nu \bar{\nu }}$ $\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ $\nu \bar{\nu }$ channel, both scaled by a factor 5000 for better visibility, and for the background processes. No selection is applied

4.2 Double Higgs production in the decay to $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$

In the $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  decay channel, the fully leptonic and semi-leptonic final states are dominated by background processes with leptons and missing transverse momentum [57]. Therefore, only the fully hadronic final state is considered here. The analysis is optimised separately for 1.4 and 3 TeV.

After the initial classification, the candidate events for $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ $\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$ are re-clustered into six jets using the longitudinally invariant $k_{t}$ algorithm with a radius parameter of $R=0.7$. The six jets are grouped by minimising

$$\begin{aligned} \chi ^{2} = \frac{(m_{ij} - m_{\text {H}})^{2}}{\sigma ^2_{\text {H}\rightarrow \mathrm{b}\bar{\mathrm{b}}}} + \frac{(m_{klmn} - m_{\text {H}})^{2}}{\sigma ^2_{\text {H}\rightarrow \text {W}\text {W}^{*}}} + \frac{(m_{kl} - m_{\text {W}})^{2}}{\sigma ^2_{\text {W}}} , \end{aligned}$$

(2)

where i, j, k, l, m, n are the indices denoting the six jets. The invariant mass resolutions are obtained from fitting a Gaussian-like function with asymmetrical width parameters to the respective peaks in the invariant mass spectra of the decay products, obtaining $\sigma _{\text {H}\rightarrow \mathrm{b}\bar{\mathrm{b}}}$ = 15.0 GeV (8.4 GeV), $\sigma _{\text {H}\rightarrow \text {W}\text {W}^{*}}$ = 36.6 GeV (7.4 GeV), and $\sigma _{\text {W}}$ = 13.1 GeV (9.5 GeV) at 3 TeV for the width below (above) the maximum, and similar values at 1.4 TeV [57]. To suppress background processes without b-quarks while minimising signal loss, the highest b-tag value among the six jets has to be at least 0.7 at 3 TeV. At 1.4 TeV, the second highest b-tag value is required to be above 0.2. As the contribution of s-channel processes such as $\text {e}^{+}\text {e}^{-}\rightarrow \mathrm{q}\bar{\mathrm{q}}\mathrm{q}\bar{\mathrm{q}}$ compared to W-boson fusion processes is larger at 1.4 TeV, in addition the total transverse momentum of the Higgs boson pair is required to be larger than 30 GeV, which enhances the fraction of processes with neutrinos in the final state [57].

The signal selection is performed using Boosted Decision Trees (BDTs) trained on the following input variables [57]: Invariant masses and angular distributions of the $\mathrm{b}\bar{\mathrm{b}}$ system, of the WW$^{*}$ system, of the jets associated with the W decay, and of the $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$ system, as well as the energy of the jets originating from the W boson are provided to the training. In addition, the transverse momenta of the two reconstructed Higgs bosons and of the di-Higgs system, angular variables between the rest-frames of the $\mathrm{b}\bar{\mathrm{b}}$, WW$^{*}$ and HH systems and of the jets associated with the W decays as well as the sphericity, the merging scales of exlusive jet clustering, and b- and c-tag values are used as input to the BDT training as well. A cut on the BDT response is applied to maximise the precision of the cross section measurement. The resulting event yields in the signal region for the HH$\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  signal and the main background processes are listed in Table 1 for 1.4 TeV and in Table 2 for 3 TeV. Although the signal selection is optimised for the decay channel $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$, there are significant contributions from other Higgs decay channels as well.

Comparing the efficiencies between the two collision energies, the Higgs bosons in the ${\text {H}\text {H}\nu \bar{\nu }}$ signal become more forward at 3 TeV. This is one of the reasons why the processes with $\text {e}^{\pm }\mathrm {\gamma }$ initial state are kinematically more similar to the signal at 3 TeV, making it more difficult to suppress them.

Table 1 Cross sections, $\sigma$, selection efficiencies, $\epsilon _{\text {BDT}}$, and expected number of events in the HH$\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  signal region, $N_{\text {BDT}}$, at $\sqrt{s} =$1.4 TeV for $\mathcal {L}=2.5$ ab$^{-1}$. The cross sections are for unpolarised beams; the number of events assumes the 4:1 polarisation scheme [16]

Table 2 Cross sections, $\sigma$, selection efficiencies, $\epsilon _{\text {BDT}}$, and expected number of events in the HH$\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  signal region, $N_{\text {BDT}}$, at $\sqrt{s} =$ 3 TeV for $\mathcal {L}=5$ ab$^{-1}$. The cross sections are for unpolarised beams; the number of events assumes the 4:1 polarisation scheme [16]

Table 3 Cross sections, $\sigma$, selection efficiencies, $\epsilon _{\text {BDT}}$, and expected number of events in the HH$\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$ signal region, $N_{\text {BDT}}$, at $\sqrt{s} =$1.4 TeV for $\mathcal {L}=2.5$ ab$^{-1}$. The cross sections are for unpolarised beams; the numbers of events assume the 4:1 polarisation scheme [16]

Table 4 Cross sections, $\sigma$, selection efficiencies, $\epsilon _{\text {looseBDT}}$ ($\epsilon _{\text {tightBDT}}$), and expected number of events in the loose (tight) BDT selection region of the HH$\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  analysis, $N_{\text {looseBDT}}$ ($N_{\text {tightBDT}}$), at $\sqrt{s} =$3 TeV for $\mathcal {L}=5$ ab$^{-1}$. The cross sections are for unpolarised beams; the numbers of events assume the 4:1 polarisation scheme

4.3 Double Higgs production in the decay to $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$

Candidate events for the final state $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  at 3 TeV are pre-selected according to the orthogonality selection (Sect. 4.1). The events are re-clustered with the VLC algorithm which provides an improvement in the di-jet mass resolution as shown in [52]. The VLC algorithm is applied in exclusive mode requiring $N=4$ jets and using a radius parameter $R=1.1$ and the parameters $\beta = \gamma = 1$. To enhance the signal fraction at $\sqrt{s} =$1.4 TeV, if $\sum _{4}{b\mathrm {-tag}} < 2.3$, events are required to have a sum of the jet energy of $\sum {E(\text {jet})}>$ 150 GeV and the second highest jet transverse momentum must be $p_{T}(\text {jet}_{2}) >$ 25 GeV.

Since both Higgs bosons are expected to be on-shell, the four jets are then grouped as two Higgs candidates by minimising the absolute difference between the resulting di-jet masses $|m_{ij} - m_{kl}|$. BDTs are trained based on the pre-selected events in order to optimise the signal selection efficiency and purity.

The following observables were chosen for the multivariate analyses: the sum of all b-tag weights, the ratio between the sum of all c-tag weights and the sum of all b-tag weights, the invariant mass of each jet pair, the cosine of the angle between the two paired jets for each jet pair evaluated in the centre-of-mass system, the total invariant mass of the system, the missing transverse momentum computed as the opposite of the vectorial sum of the momenta of all jets, the number of photons with energy larger than 25 GeV, and the maximum absolute pseudorapidity among the four jets. These ovservables are sensitive to various properties distinguishing the signal from background processes such as the presence of heavy flavour jets and neutrinos, invariant mass and angular distributions of Higgs boson decay products, and to differences of the W-boson fusion to the s-channel topology. The analyses are optimised separately for 1.4 and 3 TeV.

For the cross section measurement, the cut on the BDT response is optimised for the signal significance. The resulting expected event yields for the 1.4 TeV analysis are listed in Table 3. At 3 TeV, two selections are defined: the “tight BDT” region with a BDT cut of BDT response $\,>\,0.12$, which is optimised for signal significance, and the “loose BDT” region with a cut of BDT$\,>\,0.05$, which is optimised for the extraction of the Higgs self-coupling. The expected event yields for the two selection variants at 3 TeV for a luminosity of $\mathcal {L}=5$ ab$^{-1}$ are listed in Table 4. Both selection regions contain also a significant contribution from decays other than $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$. As the processes with $\text {e}^{\pm }\mathrm {\gamma }$ initial state do not produce more than 2 b-jets, they are strongly suppressed by criteria based on the b-tag weights. This makes the fraction of these backgrounds passing the event selection less energy-dependent than in the $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  analysis. This is also reflected in the fact that the BDT selection is more efficient for the signal in the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  analysis at both energies.

While the cross section is measured in the tight BDT region, the expected precision on ${g_{\text {H}\text {H}\text {H}}}$ and ${g_{\text {H}\text {H}\text {W}\text {W}}}$ is evaluated based on differential distributions in the loose BDT region to allow for a larger event sample. Figure 5 shows the distribution of the BDT response in the loose BDT region. From Fig. 5a, it can be seen that the SM ${\text {H}\text {H}\nu \bar{\nu }}$ signal is dominant compared to backgrounds at higher BDT score values. Selected samples with modified ${g_{\text {H}\text {H}\text {H}}}$ are compared in Fig. 5b, which shows a small overall sensitivity of the BDT score to the Higgs self-coupling. The main influence on the area of the distributions is the total cross section: The selection efficiencies vary only between 17 and 18% among the event samples with the given coupling values, while the total cross sections vary between 0.471 and 0.68 fb. The distribution of the invariant mass of the double Higgs boson system for the SM contributions in the loose BDT region is presented in Fig. 6. Figure 7a, b show the invariant di-Higgs mass distributions for selected values of ${g_{\text {H}\text {H}\text {H}}}$ and ${g_{\text {H}\text {H}\text {W}\text {W}}}$. The di-Higgs invariant mass distributions between points with similar, but opposite, variation of the ${g_{\text {H}\text {H}\text {H}}}$ coupling differ especially in the lower invariant mass region as illustrated in Fig. 7b, comparing the distributions between ${\kappa _{\text {H}\text {H}\text {H}}}= 0.8, 1.2$ and 2.2 to the SM. As shown in Fig. 7, the ${g_{\text {H}\text {H}\text {W}\text {W}}}$ coupling impacts also the higher invariant mass region, which allows it to be distinguished from modifications in the ${g_{\text {H}\text {H}\text {H}}}$ coupling.

Fig. 5

BDT response distribution of a all SM contributions stacked and b a selection of signal samples with modified ${g_{\text {H}\text {H}\text {H}}}$ in the loose BDT selection at 3 TeV CLIC. The Higgs-gauge boson coupling ${g_{\text {H}\text {H}\text {W}\text {W}}}$ is kept at its SM value

Fig. 6

Invariant mass of the Higgs boson pair for the SM contributions in the loose BDT selection at 3 TeV CLIC

Fig. 7

Invariant mass of the Higgs pairs for the signal contributions with different values of ${\kappa _{\text {H}\text {H}\text {H}}}$ and ${\kappa _{\text {H}\text {H}\text {W}\text {W}}}$ in the loose BDT region. a Comparing samples with one of the couplings fixed to the SM value. b Comparing samples with ${\kappa _{\text {H}\text {H}\text {H}}}<1$ and ${\kappa _{\text {H}\text {H}\text {H}}}>1$. The sample with ${\kappa _{\text {H}\text {H}\text {H}}}=2.2$ has roughly the same total cross section as the SM case

5 Cross section measurement

5.1 Precision of the cross section measurement for ${\text {H}\text {H}\nu \bar{\nu }}$ production at 1.4 and 3 TeV

The cross-section measurement is based on the baseline luminosity and polarisation scheme resulting in the event yields for the WBF Higgs pair production signal and the backgrounds listed in Tables 1 and 2 for the $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  analysis and in Tables 3 and 4 for the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  analysis. From this, the precision of the cross-section measurement assuming the SM value can be determined according to $\frac{\Delta \sigma }{\sigma } = \frac{\sqrt{S+B}}{S}$, where S (B) is the number of signal (background) events passing the selection. In the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  ($\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$) analysis channel, the contribution of $\text {H}\text {H}$ decaying to other final states than $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  ($\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$) yet passing the signal selection is counted towards the number of signal events.

The $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$  contribution to the ${\text {H}\text {H}\nu \bar{\nu }}$ final state exhibits a dependence on ${g_{\text {H}\text {H}\text {H}}}$, though a different one than the WBF component. In the ${\text {H}\text {H}\nu \bar{\nu }}$ analysis at 3 TeV, the modification of the ${\text {H}\text {H}\nu \bar{\nu }}$ cross section due to the variation of ${g_{\text {H}\text {H}\text {H}}}$ is treated as independent of a possible change in the $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$  contribution due to analysis selection criteria. It has been checked that the impact of a different efficiency for the $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$  component is small. In a future study, the components could be separated in the signal region based on kinematic information, using their individual dependencies on ${g_{\text {H}\text {H}\text {H}}}$.

The energy stage at $\sqrt{s}={1.4}$ TeV with an integrated luminosity of $\mathcal {L}$ = 2.5 ab$^{-1}$ and the 4:1 polarisation scheme provides evidence for the $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ process with a measurement significance of 3.5 $\sigma$ corresponding to a cross-section precision of 28%. At the 3 TeV stage alone, the observation of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ production is reached after 700 fb$^{-1}$ of data taking. Based on the 3 TeV stage and both decay channels, the precision of the ${\text {H}\text {H}\nu \bar{\nu }}$ cross-section measurement is 7.3%. The 3 TeV stage clearly dominates the cross-section measurement for WBF double Higgs production. With the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  channel at 3 TeV alone, the precision is 7.4%. This demonstrates that the contribution from the $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  analysis is very small. In the following, we therefore consider only the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  analysis. The uncertainties on the cross section measurement are summarised in Table 5.

As described in Sect. 2, the $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ cross section is dependent on the beam polarisation. In the nominal 4:1 polarisation scheme, the number of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ events is scaled by a factor of $f_p=1.43$ (1.48) at 1.4 TeV (3 TeV). For the ZHH  process at 1.4 TeV, the polarisation factor is $f_p= 1.072$. The background composition depends on the electron beam polarisation modes as well. As described in Sect. 3.1, the background kinematics have been found to be mostly independent of the polarisation. Therefore, unpolarised beams are used for the simulation and the polarisation is only taken into account in the cross section of the background processes. Some of the backgrounds scale by the same polarisation factor as the signal, others are influenced less by the polarisation. We scale all backgrounds by the same factor $f_p=1.48$. This constitutes an upper limit for the background in the negative polarisation run and a lower limit in the positive polarisation run. Since overall more luminosity is collected with negative beam polarisation, this is a conservative approach. Table 6 shows the dependence of the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  cross-section measurement uncertainty on the polarisation.

Table 5 Measurement uncertainties for the cross section of $\text {e}^{+}\text {e}^{-}\rightarrow {\text {H}\text {H}\nu \bar{\nu }}$ at the different stages of CLIC with different collision energy $\sqrt{s}$ and integrated luminosity $\mathcal {L}$ including different decay channels and assuming the 4:1 polarisation scheme

Table 6 Dependence of the cross-section measurement for ${\text {H}\text {H}\nu \bar{\nu }}$ at 3 TeV on the distribution of the luminosity between the two beam polarisation states of the electron beam. The same polarisation factor is assumed for signal and background as explained in the text

Only statistical uncertainties are considered in this study. Systematic uncertainties for the measurement of the single Higgs production cross section $\sigma ({\text {H}\nu \bar{\nu }}) \times BR(\text {H}\rightarrow \mathrm{b}\bar{\mathrm{b}})$ from various potentially dominant sources of systematic uncertainties are evaluated in [16]. Potential sources include the luminosity spectrum, the total luminosity, the beam polarisation, the jet energy scale and flavour tagging. For the $\sigma ({\text {H}\nu \bar{\nu }}) \times BR(\text {H}\rightarrow \mathrm{b}\bar{\mathrm{b}})$ measurement, they are shown to be at the per mille level. As the Higgs bosons in the ${\text {H}\text {H}\nu \bar{\nu }}\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  process are kinematically similar, the systematic uncertainties are expected to be of similar size. Compared with the almost two orders of magnitude higher statistical uncertainty, the systematic uncertainties are assumed to be irrelevant for this study.

The dependence of the cross section on the value of the trilinear Higgs self-coupling (Fig. 3) is used to derive the projected uncertainty for the extraction of the trilinear Higgs self-coupling from the measurement of the cross section.

In order to determine the expected precision for the measurement at CLIC, a template fit is used based on full detector simulation of event samples with different values of ${g_{\text {H}\text {H}\text {H}}}$ and ${g_{\text {H}\text {H}\text {W}\text {W}}}$. A $\chi ^2$ minimisation is performed, using the SM sample as the observed data. For the cases with ${g_{\text {H}\text {H}\text {W}\text {W}}}={g_{\text {H}\text {H}\text {W}\text {W}}^{\text {SM}}}$, pseudo-experiments are drawn in order to determine the confidence interval at 68% C.L. among the resulting measurements of ${g_{\text {H}\text {H}\text {H}}}$. Based on the measurement of the ${\text {H}\text {H}\nu \bar{\nu }}$ production cross section at 3 TeV only, the expected constraints at 68% C.L. for ${\kappa _{\text {H}\text {H}\text {H}}}$, assuming the SM value for ${g_{\text {H}\text {H}\text {W}\text {W}}}$, are $[0.90, 1.12] \cup [2.40,2.61]$.

5.2 Precision with ${\text {H}\text {H}\nu \bar{\nu }}$ and ${\text {Z}\text {H}\text {H}}$  production at 1.4 TeV

One approach that resolves the ambiguity on ${g_{\text {H}\text {H}\text {H}}}$ arising from the ${\text {H}\text {H}\nu \bar{\nu }}$ cross-section measurement is the combination with a measurement of the double Higgsstrahlung cross section, as described in Sect. 2. The estimates were done for $\sqrt{s}$ = 1.4 TeV as this is the energy stage of CLIC at which the ${\text {Z}\text {H}\text {H}}$  cross section is largest. No dedicated full-simulation study has been conducted. For illustration, similar analyses have been performed in full simulation for ILC at $\sqrt{s} = 500$ GeV and CLIC at $\sqrt{s} = 3$ TeV. At the ILC with $\sqrt{s} = 500$ GeV [12], a signal efficiency of 19% and a background level of 3.6 times the number of signal events is reached for the hadronic decays of the Z boson. In the CLIC study at 3 TeV [58], the signal efficiency is 20% with twice the number of background than signal events. In both cases, the numbers refer to the analyses with hadronic Z decays and the $\text {H}\text {H}\rightarrow$ $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  channel.

Based on assumptions for different signal efficiencies and background levels, Table 7 lists the significance with which the Higgsstrahlung process $\text {Z}\text {H}\text {H}$ can be observed at $\sqrt{s}$ =1.4 TeV with $\mathcal {L}$ = 2.5 ab$^{-1}$ integrated luminosity. In addition, the CLIC energy stage at 1.4 TeV will provide a cross-section measurement of the ${\text {H}\text {H}\nu \bar{\nu }}$ process in the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  final state with a luminosity of 2.5 ab$^{-1}$ and the 4:1 polarisation scheme applied, cf. Table 3. On its own, this measurement leads to the constraints [0.64, 2.3] at 68% C.L. in ${\kappa _{\text {H}\text {H}\text {H}}}$. Combining the cross-section measurements of ${\text {H}\text {H}\nu \bar{\nu }}$ and ZHH  at the CLIC stage at 1.4 TeV results in the constraints [0.71, 1.67] for a signal efficiency of 20% and a background level of twice the signal event number, which is well-motivated by the full-simulation studies described above. After the CLIC run at the 3 TeV energy stage, the differential measurement of ${\text {H}\text {H}\nu \bar{\nu }}$ production will significantly improve those constraints as described in the next section. Then, the contributions from the measurements at 1.4 TeV in ${\text {H}\text {H}\nu \bar{\nu }}$ and ${\text {Z}\text {H}\text {H}}$  will be small (cf. Table 9 and Fig. 9).

Table 7 Significance for ${\text {Z}\text {H}\text {H}}$  at 1.4 TeV in dependence on the assumptions for the performance of the ZHH analysis at 1.4 TeV. The signal efficiency $\epsilon _{\text {Sig}}$ and the ratio of selected events from background to signal, B/S, are varied

6 Self-coupling extraction based on sensitive kinematic observables

6.1 Expected precision for the trilinear Higgs self-coupling g$_{\text {HHH}}$

Differential distributions sensitive to new physics in the Higgs self-coupling can be used to measure more precisely the trilinear Higgs self-coupling ${g_{\text {H}\text {H}\text {H}}}$ and the quartic coupling to W bosons ${g_{\text {H}\text {H}\text {W}\text {W}}}$ [25]. Based on the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  selection, we make use of kinematic observables sensitive to the Higgs self-coupling as described in Sect. 2. The highest sensitivity can be reached when using the invariant mass of the Higgs boson pair in bins of the BDT score. Figure 8 shows the kinematic bins that are used for a template fit to determine the expected confidence intervals on ${g_{\text {H}\text {H}\text {H}}}$ exclusively and on ${g_{\text {H}\text {H}\text {H}}}$ and ${g_{\text {H}\text {H}\text {W}\text {W}}}$ simultaneously. The one-parameter fit in ${g_{\text {H}\text {H}\text {H}}}$ based on the differential measurement of ${\text {H}\text {H}\nu \bar{\nu }}$ at 3 TeV results in expected constraints on ${\kappa _{\text {H}\text {H}\text {H}}}$ of [0.92, 1.12] at 68% C.L.

Fig. 8

Kinematic bins used for the ${\text {H}\text {H}\nu \bar{\nu }}$ sensitivity at 3 TeV: the invariant mass of the Higgs boson pair M(HH) in bins of the BDT score

As discussed in Sect. 5.2, the influence of the ZHH measurement at the second stage is estimated using assumptions based on full simulation studies. We therefore assume in the following that a signal efficiency of 20% and a background level of twice the signal number can be achieved. This performance is applied to the full visible branching fraction of the Z boson, although leptonic decay channels have been found to give much larger signal efficiencies [12]. Several different cases of signal efficiencies and background levels are compared in Table 8 showing that the resulting uncertainties on the Higgs trilinear self-coupling are rather stable.

Figure 9 illustrates the resulting ${\Delta \chi ^{2}}$curves from the different steps of the analysis: Adding the information from the ${\text {Z}\text {H}\text {H}}$  analysis  to the rate-only measurement of ${\text {H}\text {H}\nu \bar{\nu }}$ raises the second minimum above the 68% C.L. However, with the differential measurement at 3 TeV alone, the second minimum is removed already, and the expected constraints for ${\kappa _{\text {H}\text {H}\text {H}}}$ at 68% C.L. are [0.92, 1.12], as discussed above. In this case, the impact of the ${\text {Z}\text {H}\text {H}}$  analysis at 1.4 TeV is small. By combining the differential analysis of ${\text {H}\text {H}\nu \bar{\nu }}$ with the ${\text {Z}\text {H}\text {H}}$  and ${\text {H}\text {H}\nu \bar{\nu }}$ cross-section measurements at 1.4 TeV, the best constraints are obtained, reaching [0.92, 1.11] at 68% C.L. This is the final resulting expectation for the sensitivity of the full CLIC programme to the trilinear Higgs self-coupling using the invariant di-Higgs mass and the BDT score as template. Table 9 summarises the 68% C.L. constraints obtained for ${g_{\text {H}\text {H}\text {H}}}$/${g_{\text {H}\text {H}\text {H}}^{\text {SM}}}$ with the different approaches.

Table 8 Constraints on ${\kappa _{\text {H}\text {H}\text {H}}}$ in dependence on the assumptions for the performance of the ZHH analysis at 1.4 TeV. The signal efficiency $\epsilon _{\text {Sig}}$ and the ratio of selected events from background to signal, B/S, are varied. To obtain these constraints, the ${\text {Z}\text {H}\text {H}}$  measurement is combined with the full simulation results from the ${\text {H}\text {H}\nu \bar{\nu }}$ analyses at 1.4 and 3 TeV, using differential information for ${\text {H}\text {H}\nu \bar{\nu }}$ at 3 TeV

Fig. 9

$\Delta \chi ^2$ curves based on rate-only and differential information in the ${\text {H}\text {H}\nu \bar{\nu }}$ measurement at 3 TeV without and with a combination with the measurement of the ${\text {Z}\text {H}\text {H}}$  production cross section at 1.4 TeV. As a comparison, the $\Delta \chi ^2$ for the case of the second energy stage only is shown

Table 9 Constraints on ${\kappa _{\text {H}\text {H}\text {H}}}$ obtained in the full detector simulation study using a multivariate analysis for selection. The constraint from cross section only is obtained in the tight BDT selection. The constraints based on differential distributions are derived in the loose BDT selection

These results can be interpreted in scenarios of new physics modifying the Higgs self-coupling. An example is the case of a Higgs plus singlet model [21, Sec. 6.1] and [59], where a general real singlet scalar is added to the SM Higgs sector. This could lead to a strong first-order electroweak phase transition and therefore offers an explanation of baryogenesis. The model introduces a heavy singlet as well as a mass eigenstate mixing with the SM-like Higgs boson. The new parameters are therefore the mass of the heavy scalar and the mixing angle between the singlet and the SM-like Higgs boson, as well as the parameters of the scalar potential. The parameter space is constraint by theoretical considerations such as unitarity and perturbativity. In addition, the electroweak (EW) vacuum is required to be stable. Considering direct searches in resonant double Higgs boson production as well as the sensitivity to the trilinear Higgs self-coupling, CLIC is sensitive to a sizable fraction of the parameter space which is also compatible with a strong first-order EW phase transition. In most cases, single Higgs boson coupling measurements give a similar reach as the Higgs boson pair searches, allowing a complementary assessment of the implications on the electroweak phase transition.

6.2 Expected precision for simultaneous fit of g$_{\text {HHH}}$ and g$_{\text {HHWW}}$

As described in Sect. 2, the Higgs-gauge vertex ${\text {H}\text {H}\text {W}\text {W}}$ contributes to ${\text {H}\text {H}\nu \bar{\nu }}$ as well. We can therefore extend the study of ${\text {H}\text {H}\nu \bar{\nu }}$ production at 3 TeV to fit simultaneously the modified couplings ${\kappa _{\text {H}\text {H}\text {H}}}$ and ${\kappa _{\text {H}\text {H}\text {W}\text {W}}}$. All other EFT couplings are kept at the SM value, in particular the coupling $g_{\text {Z}\text {Z}\text {H}\text {H}}$. Based on the differential distribution and binning depicted in Fig. 8, we determine the 68% and 95% C.L. contours for two degrees of freedom. The deviation of the nominal samples from the SM by $\Delta \chi ^2=2.3$ (6.18) is used for the constraints at 68% (95%) C.L.

The resulting constraints are shown in Fig. 10. At 68% C.L. the simultaneous fit leads to expected constraints of up to $20\%$ in ${\kappa _{\text {H}\text {H}\text {H}}}$ and up to 4% in ${\kappa _{\text {H}\text {H}\text {W}\text {W}}}$ across the allowed range of the other coupling. Due to the anticorrelation illustrated in Fig. 10, the individual constraints for fixed values of the other coupling are substantially smaller. Going beyond the two effective couplings ${g_{\text {H}\text {H}\text {W}\text {W}}}$ and ${g_{\text {H}\text {H}\text {H}}}$, this measurement can be combined with other measurements at CLIC in order to perform a global EFT fit of the full set of relevant operators.

Fig. 10

Confidence contours at 68% and 95% C.L. for the simultaneous fit of ${\kappa _{\text {H}\text {H}\text {H}}}$ and ${\kappa _{\text {H}\text {H}\text {W}\text {W}}}$ based on differential measurement in ${\text {H}\text {H}\nu \bar{\nu }}$ production at 3 TeV CLIC

7 Conclusions

In this paper, the prospects for the extraction of the trilinear Higgs self-coupling and the quartic ${\text {H}\text {H}\text {W}\text {W}}$ coupling at CLIC are presented. The results are based on double Higgs-boson production in the processes $\text {e}^{+}\text {e}^{-}\rightarrow$ ${\text {H}\text {H}\nu \bar{\nu }}$ and $\text {e}^{+}\text {e}^{-}\rightarrow$ ${\text {Z}\text {H}\text {H}}$. Analyses of ${\text {H}\text {H}\nu \bar{\nu }}$ production have been performed for the decay channels $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  and $\mathrm{b}\bar{\mathrm{b}}\mathrm{W}\mathrm{W}^{*}$  in full simulation. The analyses assume the second and third stage of CLIC at collision energies of 1.4 TeV and 3 TeV. In addition, the contribution of the ${\text {Z}\text {H}\text {H}}$  cross-section measurement has been included for 1.4 TeV. The channel with the highest sensitivity to the Higgs self-coupling and the ${\text {H}\text {H}\text {W}\text {W}}$ coupling at CLIC is the $\mathrm{b}\bar{\mathrm{b}}\mathrm{b}\bar{\mathrm{b}}$  decay channel of ${\text {H}\text {H}\nu \bar{\nu }}$ production at 3 TeV, where the total cross-section measurement as well as differential distributions can be used to extract the couplings. The differential measurement is based on the invariant mass distribution of the double Higgs-boson system as well as a multivariate score.

Generally, this channel can be useful to study the impact of heavy flavour tagging and jet energy resolution, especially in the forward direction, realised in the CLIC detector models. In this case, the CLIC_ILD model was used. No significant change is expected for the application of this analysis to the current CLICdet model. This analysis benefits from the higher centre-of-mass energy due to the increase in cross section of ${\text {H}\text {H}\nu \bar{\nu }}$ production. It therefore provides a strong motivation for the CLIC 3 TeV energy stage.

Beams in the simulated samples used in this analysis are unpolarised. Based on the cross section for each polarisation mode, the results are scaled to the baseline polarisation scheme of CLIC with the running time shared between P(e$^{-}$) = $-80\%~(+\,80\%)$ in the ratio 4:1. Future studies making use of the polarisation-dependent kinematic behavior might improve the signal selection. Furthermore, future studies could treat the $\text {Z}(\rightarrow \nu \bar{\nu })\text {H}\text {H}$  component separately in the ${\text {H}\text {H}\nu \bar{\nu }}$ final state. This contribution is particularly important at 1.4 TeV and also depends on the polarisation mode. As it also has a dependence on ${g_{\text {H}\text {H}\text {H}}}$, albeit a different one than the WBF production channel, this can be exploited separately.

At the 1.4 TeV energy stage of CLIC, evidence for the ${\text {H}\text {H}\nu \bar{\nu }}$ process of the SM can be reached with a significance of 3.5 $\upsigma$. With a luminosity of only 700 fb$^{-1}$, the process can be observed with 5.0 $\upsigma$ at 3 TeV. Taking into account only the 1.4 TeV stage of CLIC with cross-section measurements of ${\text {H}\text {H}\nu \bar{\nu }}$ and ZHH  allows the measurement of the Higgs self-coupling ${g_{\text {H}\text {H}\text {H}}}$ with relative uncertainties of $-29\%$ and $+67\%$ around the SM value at 68% C.L. Based on events of double Higgs-boson production at both high-energy stages, CLIC can be expected to measure the trilinear Higgs self-coupling ${g_{\text {H}\text {H}\text {H}}}$ with a relative uncertainty of $-8\%$ and $+11\%$ at 68% C.L., assuming the Standard Model and setting the quartic ${\text {H}\text {H}\text {W}\text {W}}$ coupling to its Standard Model value. Measuring simultaneously the trilinear Higgs self-coupling and the quartic Higgs-gauge coupling results in constraints at 68% C.L. below 4% in ${g_{\text {H}\text {H}\text {W}\text {W}}}$ and below 20% in ${g_{\text {H}\text {H}\text {H}}}$ for large modifications of ${g_{\text {H}\text {H}\text {W}\text {W}}}$. These results illustrate the strength of the proposed CLIC programme to make a precise measurement of the trilinear Higgs self-coupling.