Search for the $HH \rightarrow b \bar{b} b \bar{b}$ process via vector-boson fusion production using proton-proton collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search for Higgs boson pair production via vector-boson fusion (VBF) in the $b\bar{b}b\bar{b}$ final state is carried out with the ATLAS experiment using 126 fb$^{-1}$ of proton-proton collision data delivered at $\sqrt{s} = 13$ TeV by the Large Hadron Collider. This search is sensitive to VBF production of additional heavy bosons that may decay into Higgs boson pairs, and in a non-resonant topology it can constrain the quartic coupling between the Higgs bosons and vector bosons. No significant excess relative to the Standard Model expectation is observed, and limits on the production cross-section are set at the 95 % confidence level for a heavy scalar resonance in the context of an extended Higgs sector, and for non-resonant Higgs boson pair production. Interpretation in terms of the coupling between a Higgs boson pair and two vector bosons is also provided: coupling values normalised to the Standard Model expectation of $\kappa_{2V}<-0.43$ and $\kappa_{2V}>2.56$ are excluded at the 95 % confidence level in data.


Introduction
The Higgs boson ( ) was discovered by the ATLAS and CMS collaborations in 2012 [1, 2] using proton-proton ( ) collisions at the Large Hadron Collider (LHC). The measured properties have so far been found to be in agreement with the Standard Model (SM) predictions. The production of a pair of Higgs bosons ( ) is a rare process in the SM with a cross-section about 1000 times smaller than the single Higgs boson production cross-section, but various theories beyond the SM (BSM) predict cross-sections for production that are significantly higher than the SM prediction. Spin-0 resonances, with narrow or broad width, that decay into Higgs boson pairs, appear in BSM scenarios [3,4]. Enhanced non-resonant Higgs boson pair production is predicted by many models, for example those featuring light coloured scalars [5] or new contact interactions, such as direct¯vertices [6,7].
Previous searches for Higgs boson pair production in the¯¯channel were carried out in the gluon-gluon fusion (ggF) production mode by the ATLAS and CMS collaborations [8][9][10][11][12], and limits were set for resonant and non-resonant production. Statistical combinations of search results for in various decay channels were also performed by the two experiments [13,14], profiting from the sensitivity of several final states. This paper focuses on searches for Higgs boson pair production via vector-boson fusion (VBF), through diagrams such as those presented in Figure 1, and using the dominant →¯decay mode [15]. The VBF process ( → ) is characterised by the presence of two jets ( ) with a large rapidity gap resulting from quarks from which a vector boson ( ) is radiated. In the SM, three different types of couplings are involved in production via VBF: the Higgs boson self-coupling ( ), the Higgs-boson-vector-boson coupling ( ) and the quartic (di-vector-boson-di-Higgs-boson, or ) coupling. The coupling modifiers , and 2 control the strength of the , and couplings with respect to the SM value, respectively, and are normalised so that they are equal 1 in the SM. A deviation of these coupling modifiers from their SM expectations could lead to enhanced production. While searches in the ggF mode are more sensitive to deviations in , the VBF topology has unique sensitivity to 2 [15] because the ggF mode does not involve the interaction. For resonant production, two classes of signals are tested to perform a generic inclusive search for resonances with masses in the range 260-1000 GeV. The first signal class is representative of a broad resonance with width typically 10-20% of the signal mass; it corresponds to a heavy scalar of the 2HDM Type II model [16] and is obtained by setting the ratio of vacuum expectation values of the two Higgs doublets tan( ) = 2.0 and sin( − ) = 0.6, where is the mixing angle between the two CP-even Higgs bosons. The second class features a narrow resonance with a fixed width of 4 MeV.

ATLAS detector
The ATLAS experiment [17][18][19] at the LHC operates a multipurpose particle detector with a forwardbackward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadronic calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range | | < 2.5. It consists of silicon pixel, silicon microstrip, and transition-radiation tracking detectors. Lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. A hadronic steel/scintillator-tile calorimeter covers the central pseudorapidity range (| | < 1.7). The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to | | = 4.9. The muon spectrometer surrounds the calorimeters and incorporates three large air-core toroidal superconducting magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 Tm across most of the detector. The muon spectrometer includes a system of precision tracking chambers and fast detectors for triggering [20]. A two-level trigger system is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information to reduce the accepted rate to at most 100 kHz. This is followed by the software-based high-level trigger (HLT), which reduces the accepted event rate to 1 kHz on average.

Data and simulated samples
This search is performed using data collected by the ATLAS experiment between 2016 and 2018 in √ = 13 TeV LHC collisions, which correspond to an integrated luminosity of 126 fb −1 . Only events recorded during stable beam conditions and when the detector components relevant to the analysis were operating normally are considered [21]. During the 2016 data-taking, a fraction of the data was affected by an inefficiency in the vertex reconstruction in the HLT, which reduced the efficiency of the algorithms used to identify jets originating from -hadron decays; those events were not retained for further analysis. This reduces the integrated luminosity of the 2016 dataset to 24.3 fb −1 .
Simulated Monte Carlo (MC) event samples are used to model signal production and the backgrounds from top-quark-pair (¯) production. The dominant process arises from multĳet production, and is modelled using data-driven techniques.
Events with a generic scalar resonance produced via VBF and decaying into →¯¯were generated with P -B v2 [22][23][24] interfaced to P 8.186 [25] for parton showering and hadronisation, with the Higgs boson mass fixed to 125 GeV [26]. The NNPDF23LO parton distribution function (PDF) 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of set [27] and the A14 set of tuned parameters [28] for underlying-event simulation were used. Resonant signal samples were produced with masses ranging from 260 GeV to 1000 GeV.
Non-resonant production of Higgs boson pairs was simulated with M G 5_ MC@NLO [29]. The Higgs boson self-coupling and the couplings of the Higgs boson to vector bosons were set to their SM values, while the 2 coupling modifier was varied. Samples with 2 equal to 0.0, 0.5, 1.0, 1.5, 2.0, and 4.0 were generated at leading order (LO). Interference between various diagrams contributing to the non-resonant signal is considered in the simulation. A linear combination of three samples is used to derive distributions for a finer granularity of 2 values, following a technique used to generate distributions [13]. A particular choice of three samples is done based on the 2 range of the generated distributions to avoid large weights and to reduce statistical uncertainties. The cross-section of the VBF process, evaluated at next-to-next-to-next-to-leading order (N 3 LO) in QCD, is 1.73 ± 0.04 fb in the SM [30][31][32][33]. The N 3 LO to LO cross-section ratio at the SM value is calculated and this factor is applied to the cross-sections at each 2 point.
To estimate the contribution from ggF production with two additional jets that can mimic the VBF topology, the SM non-resonant production of Higgs boson pairs via ggF was simulated with M G 5_ MC@NLO using the CT10 PDF set [34] and the FTapprox method [35] to include finite top-quark mass effects. In this sample, the generation of → + parton is done at next-to-leading order (NLO). Parton showers and hadronisation were simulated with Herwig 7.0.4 [36]. Interference effects with other SM processes are found to be marginal and are ignored. The cross-section is evaluated at next-to-next-to-leading order (NNLO) with the resummation at next-to-next-to-leading-logarithm (NNLL) accuracy and including top-quark mass effects at NLO [35,[37][38][39][40][41][42]; it is equal to 31.05 +1. 40 −1.99 fb. The uncertainty includes the variations of the factorisation and renormalisation scales, PDF and S .
The generation of¯events was performed with P -B v2 [43] using the NNPDF3.0NLO [44] PDF set. The parton showers, hadronisation, and underlying event were simulated using P 8.230 [45] with the NNPDF23LO PDF set and the corresponding A14 set of tuned underlying-event parameters. The predicted¯production cross-section is 831.8 +19.8 −29.2 ± 35.1 pb as calculated with the Top++ 2.0 program to NNLO in perturbative QCD, including soft-gluon resummation to NNLL accuracy [46], and assuming a top-quark mass of 172.5 GeV. The first uncertainty comes from the independent variations of the factorisation and renormalisation scales, while the second one is associated with variations in the PDF and S , following the PDF4LHC prescription with the MSTW2008 68% CL NNLO, CT10 NNLO and NNPDF2.3 5f FFN PDF sets [27,[47][48][49].
For all simulated events, -hadron and -hadron decays were handled by E G 1.2.0 [50]. To simulate the impact of multiple interactions that occur within the same or nearby bunch crossings (pile-up), minimum-bias events generated with P 8.186 using the NNPDF2.3LO set of PDFs and the A3 set of tuned parameters [51] were overlaid on the hard-scatter process. The detector response was simulated with G 4 [52,53], and the events were processed with the same reconstruction software as that was used for the data.
Jets are reconstructed from three-dimensional topological clusters of energy deposits in the calorimeter [54] with the anti-algorithm [55] implemented in the F J package [56] with radius parameter = 0.4. Clusters are calibrated at the EM scale [57] and their energy is corrected for additional energy deposits from pile-up interactions using an area-based correction [58]. Subsequently, calibration using T -and -dependent factors derived from simulation is applied, followed by the global sequential calibration [57]. The latter reduces the flavour dependence of the calibration and energy leakage effects. The final calibration is based on in situ measurements in collision data [57]. To preferentially reject jets originating from pile-up interactions, a multivariate classification algorithm (jet vertex tagger) based on tracking information [59] is used for jets with T < 60 GeV and | | < 2.4. The selected working point provides an inclusive hard-scatter process efficiency of about 97% in that kinematic region. The efficiency in the simulation is corrected to match that measured in data. Events having jets consistent with noise in the calorimeter or non-collision backgrounds are vetoed [57].
Jets containing -hadrons are identified using a multivariate algorithm (MV2c10) [60,61], which exploits information about the jet kinematics, the impact parameters of tracks associated with the jet and the presence of displaced vertices to form the decision. The -tagging requirements result in an efficiency of 70% for jets with T > 20 GeV containing -hadrons, and the misidentification rate is 0.3% (11.2%) for light-flavour (charm) jets. These were determined in a sample of simulated¯events. For all simulated events the -tagging efficiencies are corrected to match those measured in data [60,62,63].
To further correct the -jet energy for effects that are not considered in the default calibration, a jet energy regression is used. The method uses a boosted decision tree (BDT) algorithm implemented in TMVA [64]. The BDT training is performed using variables the -jet energy resolution is sensitive to: MV2c10 score, energy leakage outside the jet cone, pile-up contamination, hard radiation from the original parton, and energy loss through semileptonic -hadron decays. Both the training and the validation are performed with simulated¯samples, resulting in approximately 10% improvement in the jet energy resolution. The performance of the jet energy regression is validated in (→ ) + /¯events in data and no mismodelling is found. The →¯mass peak is found to be closer to 125 GeV and the standard deviation divided by the mean of the mass distribution is improved by about 25% for the = 600 GeV signal sample.

Event selection
Events are selected using a combination of -jet triggers, with the lowest jet transverse energy, T , threshold at 35 GeV, jet | | < 2.5 and one or two -tagged jets. The -jet trigger efficiency is measured in data, and the simulated events are corrected to match the measured trigger efficiency.
To select events compatible with VBF production of Higgs boson pairs decaying into four -quarks, exactly four central -tagged jets with T > 40 GeV and | | < 2.0 and at least two forward jets with T > 30 GeV and | | > 2.0 are required. Events with more than four central jets are vetoed in case more than four jets fulfil the -tagging requirement. To form Higgs boson candidates, -tagged jets are used, and the forward jets are considered as possible remnants of the VBF process. The Higgs boson reconstruction procedure is the same as the one described in Ref.
[8], except that the usage of the jet energy regression changes the numerical values in the signal region definition described below. A summary of the selection criteria is provided in Table 1. Table 1: Summary of the selection criteria for capturing the VBF topology, identifying →¯¯decays, and suppressing background events. Possible remnants of the VBF process are identified using the two highest-T forward jets. Labels "lead" and "subl" refer to the leading and subleading Higgs boson candidates (ordered in T ), respectively. The definitions of the different analysis regions are also provided. The transverse momenta and masses are expressed in GeV.

VBF-jets selection
The two highest-T forward jets with opposite sign of are considered as remnants of the VBF production process if the absolute value of the pseudorapidity separation between them, Δ VBF , exceeds 5.0 and their invariant mass, VBF , is greater than 1000 GeV.

Signal kinematics selection
The four central -tagged jets are considered in three possible combinations of two-jet pairings. Their invariant mass, 4 , is used to define criteria to select signal-like events. The following Δ requirements must be satisfied; these reflect the correlation of 4 with the Lorentz boost of the Higgs bosons and the angle between their decay products in the laboratory frame: where Δ lead and Δ subl are the angular distances between the jets that form, respectively, the leading and subleading Higgs boson candidates (ordered in T ). These criteria are optimised for both non-resonant and resonant Higgs boson pair production, and the numerical values are chosen to maximise the signal sensitivity.
Out of the possible pairings fulfilling the previous selection, the combination that leads to pairs with a dĳet mass closest to the SM Higgs boson mass should be the optimal choice. However, due to energy loss through semileptonic -hadron decays, this criterion is relaxed. The mass values of 123.7 GeV for the leading Higgs boson candidate and 116.5 GeV for the subleading Higgs boson candidate are found to maximise the signal significance for a resonance with a mass of 600 GeV, which lies in the middle of the covered mass range. The same target values are used for all other signal hypotheses. For a given pairing, the quantity that corresponds to the distance of the leading and subleading Higgs boson candidate masses, in the ( lead 2 , subl 2 ) plane, from the line connecting (0 GeV, 0 GeV) and (123.7 GeV, 116.5 GeV), can be computed as: The pairing with the smallest value of is chosen. Studies based on simulation indicate that for SM non-resonant production the correct pairs are identified in at least 83% of the signal events, while for broad resonances the corresponding fraction is greater than 91%.

Selection for background suppression
In order to enhance the sensitivity to signal, various requirements are applied to suppress the background. The magnitude of the vector sum of the transverse momenta of the selected four -jets and the two VBF-jets tends to peak at lower values for signal events than for multĳet events. Consequently, it is required to be less than 60 GeV. The pseudorapidity difference between the reconstructed Higgs boson candidates, |Δ |, is required to be below 1.5, and mass-dependent requirements on the transverse momenta of the leading and subleading Higgs boson candidates, respectively lead T, and subl T, , are: where 11.6 GeV and 18.1 GeV are the experimental widths of the simulated leading and subleading Higgs boson candidates, respectively. The mass resolution of the subleading Higgs boson candidate is worse because it is composed of the lower T jets. These values are derived using a 600 GeV resonant signal sample and are similar for other signal samples.
Additional requirements are imposed to reduce the number of hadronically decaying¯events by vetoing candidate events compatible with a top-quark decay. The jet with the highest -tagging score is considered as the -jet originating from a top-quark decay and the remaining central jets are considered to stem from the -boson decay. Since the top quarks are expected to be produced centrally, only central jets are tested. All possible two-jet combinations in the event are tested and the selected combination is the one with the smallest value of , defined as: where and are the reconstructed invariant masses of the -boson and top-quark candidates, respectively. The event is vetoed if ≤ 1.5. This requirement reduces the¯contamination by about 50% with negligible impact on the signal efficiency.
All requirements listed above define the signal region (SR). The number of selected signal events divided by the number of generated events after each selection step (cumulative acceptance times efficiency) is shown in Figure 2 for the non-resonant signal as a function of the 2 coupling modifier and for the resonant signal models as a function of the generated mass. The acceptance times efficiency increases as a function of the resonance mass, while for the non-resonant signal a significant drop is observed at 2 = 1. The trigger and jet selection requirements cause the drop for 2 values around 1, while the smaller acceptance times efficiency for low-mass resonances is due to the softer T spectrum of -jets.
To estimate the background and to validate the background estimation technique, two regions orthogonal to the SR are used: the sideband region (SB) and the validation region (VR). The events in the SB and VR

Background estimation
After the event selection described in Section 5, the background is dominated by multĳet and¯events. The multĳet events constitute about 95% of the total background and are modelled using data. The remaining 5% are¯events, which are modelled using simulation. The normalisation of the all-hadronic¯background is determined from data, whereas the non-all-hadronic¯background is normalised to the SM prediction. In the SM, the contribution of the pairs produced via ggF is small compared to other backgrounds and three times larger than for the VBF production. Thus production via ggF is treated as a background in this analysis and is fixed to the SM prediction. Other minor backgrounds with contributions below 0.5% are neglected. The background estimation technique is the same as in Ref. [8].

Multĳet background
The data-driven multĳet background estimation uses data events with lower -jet multiplicity and reweights them to model events with higher -jet multiplicity. The multĳet events are selected using the same trigger and selection requirements as those used in the SR, except for the -tagging requirement. In particular, the SR requires at least four -jets ("four-tag sample"). To derive a background estimate for this region, events with at least four central jets, but with only two of them -tagged ("two-tag sample"), are used.
The events in the two-tag sample are reweighted by applying a product of two event weights. The first event weight corrects for the additional -tagged jet activity and the second event weight corrects for the kinematic differences caused by requiring additional -tagged jets. These differences can arise for a variety of reasons: the -tagging efficiency varies as a function of jet T and ; the various multĳet processes contribute with different fractions in each sample; and the fraction of events accepted by each trigger path changes. The reweighting is performed using one-dimensional distributions and is iterated until the weights converge to stable values. Details of the reweighting procedure can be found in Ref.
[8]. The weights are derived in the SB using the procedure described above and validated in the VR.

6.2¯background
The shape of the¯background is modelled using simulation. The¯events are expected to contain two -jets from the decay of two top quarks and additional jets stemming from the hadronic -boson decay or additional quarks or gluons produced together with two top quarks. To reduce the statistical uncertainty, simulated¯events in the two-tag region, corrected by the kinematic weights derived for the multĳet background, are also used. The procedure is validated in the SB, and good agreement is observed between the corrected two-tag sample and the four-tag sample within the uncertainties. Samples for all-hadronic and non-all-hadronic¯decays are handled separately.

Background normalisation
The normalisations of multĳet and all-hadronic¯backgrounds are derived simultaneously by fitting the distribution to data in the SB. The distribution differs for the two backgrounds: the region of < 0.75 is enriched in all-hadronic¯events and the region > 0.75 is enriched in multĳet events. The normalisation of the non-all-hadronic¯background is fixed to its SM prediction in the fit due to its small contribution to the total yields. Two parameters are used in the normalisation fit: multĳet and all-had.¯. The multĳet parameter scales the multĳet yield from the two-tag to the four-tag sideband region after the reweighting described in Section 6.1. The all-had.¯p arameter corrects the normalisation of all-hadronic¯yields in the four-tag sideband region. The fit results are cross-checked in the VR, where the same results are obtained within statistical uncertainties.

Systematic uncertainties
Background normalisation uncertainties are propagated from the fit, described in the previous section, which determines the multĳet and all-hadronic¯yields. The statistical uncertainty of the multĳet and all-hadronic¯normalisation parameters is accounted for, including their correlations. Two nuisance parameters are defined in the final fit described in Section 8 by calculating two eigenvectors from the covariance matrix of the normalisation fit. Furthermore, the normalisation of the multĳet background estimate is verified in the VR, where agreement with data is found; its statistical uncertainty in the VR is thus applied as a normalisation systematic uncertainty of the multĳet background.
Two shape uncertainties for the multĳet background modelling are evaluated using data from the VR. The multĳet modelling uncertainty is related to the level of agreement between the data and the background model in this region. To evaluate this uncertainty, the 4 distribution is split into low and high mass regions at 400 GeV and a linear fit to the ratio of this distribution between data and the sum of all backgrounds is performed in both mass regions. Studies show that the ratio, which appears to be constant, can be sufficiently well described by a straight-line fit, and that the change in the result is marginal for other splitting points. The +1 variation of the slope of the fitted line and the +1 variation of the fitted line with inverted slope are used as up and down variations of the uncertainty in the background shape. The yield of the varied multĳet templates is fixed to its nominal value. The varied templates are derived separately for the low-and high-4 mass regions, the corresponding systematics are treated as uncorrelated in the final fit. The kinematic reweighting uncertainty is assessed by deriving an alternative multĳet template using the same procedure as in the nominal case, but using data from the VR. This difference between the multĳet templates derived in the SB and the VR is symmetrised around the nominal template, keeping the yield fixed to its nominal value.
The shape modelling uncertainty of the¯background is evaluated by comparing the nominal¯template derived using the two-tag sample as described in Section 6 and a¯template derived using the four-tag sample after the basic preselection requirements defined in Section 5. A straight-line fit to the ratio of the two-tag to four-tag 4 templates is performed. With a procedure similar to the one applied for the multĳet shape uncertainties, the up and down variations around the nominal¯template are extracted using a straight-line fit to the +1 variation of the slope of the fitted line and the +1 variation of the fitted line with inverted slope, while keeping the yield fixed to its nominal value. This uncertainty is derived separately for the non-all-hadronic (using MC) and all-hadronic (using data)¯samples and is not correlated in the fit between these two samples.
Theoretical uncertainties in the ggF background yield are evaluated by varying the renormalisation and factorisation scales and from the uncertainty associated with the choice of PDF set. The resulting variation of the expected ggF background yield is about 10%. When considering the same sources of theoretical uncertainty for the VBF signal, its acceptance times efficiency varies by 3%.
The experimental uncertainties listed below affect only MC samples. The uncertainties in the jet energy resolution and scale are evaluated at √ = 13 TeV using in situ measurement techniques described in Ref. [57]. The sources of uncertainty in these measurements are treated as fully correlated between the T and mass scales. The resolution uncertainty is evaluated in measurements documented in Ref. [65] and is assessed by applying an additional smearing to these observables. The flavour tagging efficiency and its uncertainty for -and -jets is estimated in¯events, while the light-jet misidentification rate and uncertainty is determined using dĳet events [60,62,63]. In addition, an uncertainty in the -jet trigger efficiency is derived from the per-jet online -tagging measurements [66]. The uncertainty in the integrated luminosity is 1.7% [67], obtained using the LUCID-2 detector [68] for the primary luminosity measurements.

Results
Following the statistical procedures outlined in Ref.
[1], a test statistic based on the profile likelihood ratio [69] is used to test hypothesised values VBF of the cross-section of the signal model in units of fb. This test statistic extracts the information about the signal cross-section from a likelihood fit to the data. The likelihood function includes all parameters which describe the systematic uncertainties and their correlations discussed in Section 7. As no significant excess over the background prediction is observed, exclusion limits are computed using the asymptotic formula [69]. The exclusion limits are based on the CL s method [70], where a value of VBF is regarded as excluded at the 95% confidence level (CL) when CL s is smaller than 5%. The accuracy of the asymptotic approximation is verified with sampling distributions generated using pseudo-experiments.
The mass of the four selected -jets, 4 , in the SR is used as the final discriminant for limit setting. Figure 4 shows the distribution of data and the SM background after the background-only fit in the SR and VR. In addition, the signal prediction for the narrow-width resonance hypothesis with = 800 GeV and the non-resonant signal at 2 = 3 are shown in the SR.
Upper limits on the cross-section are set for all tested models. Figure 5 shows the 95% CL upper limits for resonant production via VBF as a function of the resonance mass for narrow-and broad-width resonance hypotheses. The significance of the excess over the background-only prediction is quantified using the local 0 -value, defined as the probability of the background-only model to produce a signal-like fluctuation at least as large as that observed in the data. The most extreme 0 -value corresponds to a local significance of 1.5 standard deviations at 550 GeV.   Figure 6: Observed and expected 95% CL upper limits on the production cross-section for non-resonant production via VBF as a function of the di-vector-boson-di-Higgs-boson coupling modifier 2 . The theory prediction of the cross-section as a function of 2 is also shown. More details on the predicted cross-section can be found in Section 3.

ATLAS
The expected and observed limits on SM non-resonant production via VBF are given in Table 2. Limits are also calculated as a function of 2 , as presented in Figure 6. The observed excluded region corresponds to 2 < −0.43 and 2 > 2.56, while the expected exclusion is 2 < −0.55 and 2 > 2.72. For 2 values deviating from the SM prediction, growing non-cancellation effects result in a harder spectrum, and thereby higher-T -jets, which in turn lead to increased signal acceptance times efficiency as shown in Figure 2. This search is therefore not sensitive to the region close to the SM prediction, corresponding to 2 = 1 .  Table 3 summarises the relative impact of the uncertainties on the best-fit signal cross-section for two different narrow-width resonance production hypotheses, with masses equal to 300 GeV and 800 GeV. Only major sources of systematic uncertainty are quoted along with the impact of the statistical uncertainty. The uncertainties of similar nature are grouped into unique categories and the fit is performed independently Table 3: Dominant relative uncertainties in the best-fit signal cross-section best fit VBF ( → → ) of hypothesised resonant signal production. The leading sources of systematic uncertainty, the total systematic uncertainty and the data statistical uncertainty are provided. Two mass points are selected: = 300 GeV with the best-fit cross-section of 140 fb and = 800 GeV with 4.7 fb, which correspond to the low and high mass regions. The groups of uncertainties do not add up in quadrature to the total uncertainty, because only the dominant uncertainties are shown and also due to correlations between the uncertainties.

Conclusion
A search for both resonant and non-resonant production of pairs of Standard Model Higgs bosons via vector-boson fusion has been carried out in the¯¯channel. The analysed data were collected from √ = 13 TeV proton-proton collisions by the ATLAS detector at the LHC in 2016-2018 and correspond to an integrated luminosity of 126 fb −1 . Results for resonant production are reported in the mass range 260-1000 GeV. The largest deviation from the background-only hypothesis is observed at 550 GeV with a local significance of 1.5 standard deviations. Upper limits on the production cross-section are set for narrow-and broad-width scalar resonances at 95% CL. Limits are also set on the cross-section for non-resonant production, and as a function of the di-vector-boson-di-Higgs-boson coupling modifier, 2 . The observed 95% CL upper limit on the SM non-resonant production cross-section is 1460 fb, compatible with the expected limit at a level below two standard deviations. The observed excluded region corresponds to 2 < −0.43 and 2 > 2.56, while the expected exclusion is 2 < −0.55 and 2 > 2.72.