Probing the Higgs self-coupling through double Higgs production in vector boson scattering at the LHC

In this work we explore the sensitivity to the Higgs self-coupling $\lambda$ in the production of two Higgs bosons via vector boson scattering at the LHC. Although these production channels, concretely $W^+W^- \to HH$ and $ ZZ \to HH$, have lower rates than gluon-gluon fusion, they benefit from being tree level processes, being independent of top physics and having very distinctive kinematics that allow to obtain very clean experimental signatures. This makes them competitive channels concerning the sensitivity to the Higgs self-coupling. In order to give predictions for the sensitivity to this coupling, we first study the role of $\lambda$ at the subprocess level, both in and beyond the Standard Model, to move afterwards to the LHC scenario. We characterize the $pp\to HHjj$ case first and then provide quantitative results for the values of $\lambda$ that can be probed at the LHC in vector boson scattering processes after considering the Higgs boson decays. We focus mainly in $pp\to b\bar{b}b\bar{b}jj$, since it has the largest signal rates, and also comment on the potential of other channels, such as $pp\to b\bar{b}\gamma\gamma jj$, as they lead to cleaner, although smaller, signals. Our whole study is performed for a center of mass energy of $\sqrt{s}=14$ TeV and for various future expected LHC luminosities.


Introduction
The observation of the Higgs boson by the ATLAS and CMS experiments [1,2] in 2012 confirmed the prediction of the last particle of the Standard Model (SM) of fundamental interactions. Although this discovery allowed us to answer many important and well established questions about elementary particle physics, it also posed a lot of new mysteries concerning the scalar sector of the SM.
One of these mysteries is that of the true value of the Higgs self-coupling λ, involved in trilinear and quartic Higgs self-interactions, and its relation to other parameters of the SM. Particularly, understanding and testing experimentally the relation between λ and the Higgs boson mass, m H , will provide an excellent insight into the real nature of the Higgs particle. This relation, given in the SM at the tree level by m 2 H = 2v 2 λ, with v = 246 GeV, arises from the Brout-Englert-Higgs (BEH) mechanism [3][4][5][6], so to really test this theoretical framework one needs to measure λ independently of the Higgs mass. Unfortunately, the value of the Higgs selfcoupling has not been established yet with precision at colliders, but there is (and will be in the future) a very intense experimental program focused on the realization of this measurement (for a review, see for instance [7][8][9][10][11]).
The Higgs trilinear coupling can be probed in double Higgs production processes at the LHC, process that have been extensively studied both theoretically in , and experimentally in [49][50][51][52][53][54]. At hadron colliders, these processes can take place through a variety of production channels, being gluon-gluon fusion (GGF) and vector boson scattering (VBS), also called vector boson fusion (VBF) in the literature, the main ones regarding the sensitivity to the Higgs selfcoupling. Focusing on the LHC case, on which we will base our posterior study, the dominant contribution to double Higgs production comes from GGF, which for √ s = 14 TeV is about a factor 17 larger than from VBS [26]. Because of this, most of the works present nowadays in the literature focus on this particular H H production channel, GGF, to study the sensitivity to λ. In fact, all these works and the best present measurement at the LHC have made possible to constraint this parameter in the range λ ∈ [−8.2, 13.2] · λ SM at the 95% CL [54].
Although GGF benefits from the highest statistics and rates, it suffers the inconveniences of having large uncertainties, being a one loop process initiated by gluons, and being dependent of the top Yukawa coupling. Double Higgs production via VBS [8,11,18,21,26,27,30,32,40] is, in contrast, a tree level process not initiated by gluons and it is independent of top physics features, leading therefore to smaller uncertainties. Also, at a fundamental level, it is well known that VBS processes involving longitudinally polarized gauge bosons, like the process V L V L → H H that we are interested in, probe genuinely the self interactions of the scalar sector of the SM. This would happen specially at high energies, such as those available at the LHC, since, in this regime, each V L behaves as its corresponding would-be-Goldstone boson φ. Therefore, testing V L V L → H H is closely related to testing φφH H interactions. In this way, a new window, qualitatively different than GGF, would be open with VBS to test λ, meaning that being able to measure these processes for the first time will be a formidable test of the SM itself, and it could even lead to the discovery of physics beyond the Standard Model (BSM). Moreover, the VBS production channel is the second largest contribution to Higgs pair production, and the VBS topologies have very characteristic kinematics, which allow us to select these processes very efficiently as well as to reject undesired backgrounds. In fact, the selection techniques for VBS configurations at the LHC have experienced a great development in the last years and have improved considerably, especially in the context of electroweak (EW) vector boson scattering V V → V V [11,[55][56][57][58]. Thus, in summary, VBS double Higgs production might be very relevant to study the sensitivity to the Higgs self-coupling, despite the fact that it is considerably smaller in size than GGF, since it could lead to a cleaner experimental signal. Besides, it will be a complementary measurement to that of GGF and will, in any case, help to improve the determination of this λ coupling with better precision.
In this work, motivated by the above commented advantages, we analyze in full detail Higgs pair production via VBS at the LHC to probe the Higgs self-coupling. To this end, we first explore and characterize the subprocesses of our interest, V V → H H with V = W, Z, both for the SM with λ = λ SM and for BSM scenarios with λ = κ λ SM , and consider values of κ between 10 and −10. For this study, we fix m H to its experimental value, m H = 125.18 ± 0.16 GeV [59], and set the Higgs vacuum expectation value (vev) to v = 246 GeV. In this way, studying the sensitivity to λ in VBS will provide the desired independent test of this coupling.
Once we have deeply studied double Higgs production at the subprocess level, we then explore in this work the LHC scenario. First we analyze the process pp → H Hjj , to fully understand the properties of this scattering, and then we study and give quantitative results for the sensitivity to the Higgs self-coupling after the Higgs decays. The production of H Hjj at the LHC, including VBS and GGF, has been studied previously in [24,33], where they focus on bbττ jj final states. Our main study is performed, in contrast, in the four bottoms and two jets final state, pp → bbbbjj , since it benefits from the highest rates. We also make predictions for the interesting pp → bbγ γjj process which, although with lower rates, leads to cleaner signatures. We would like to point out that all computations and simulations are performed at the parton level with no hadronization or detector response simulation taken into account, since the work is aimed to be a first and simple approximation to the sensitivity to λ in VBS processes at the LHC.
The paper is organized as follows: In Section 2 we study VBS double Higgs production at the subprocess level in and beyond the SM. Afterwards, in Section 3, we move on to the LHC case, exploring first the pp → H Hjj scattering in Subsection 3.1 and considering later the Higgs decays, both leading H → bb decay and subleading H → γ γ one. Then, we study both signal and background rates for pp → bbbbjj in Subsection 3.2 and pp → bbγ γjj in Subsection 3.3, providing our results for the sensitivity to λ in VBS Higgs pair production at the LHC for a center of mass energy of √ s = 14 TeV and for different and future expected luminosities. Section 4 summarizes our main conclusions.

Double Higgs production in vector boson scattering
As already stated in the paragraphs above, we are interested in exploring the sensitivity to the Higgs self-coupling, λ, through VBS processes, in particular at the LHC. For that purpose, we have to study and characterize first the subprocess that leads to the specific signal we will be dealing with once we perform the full collider analysis. This subprocess will be, in our case, the production of two Higgs bosons in the final state from the scattering of two EW gauge bosons, V V → H H , with V = W, Z. Within this context, in this section we aim to understand the role of the Higgs trilinear coupling in the SM and beyond, as well as the generic characteristics of the scattering processes W + W − → H H and ZZ → H H .
The Higgs self-coupling is only present, at the tree level and in the Unitary gauge, in the s-channel diagram of the studied processes, so the sensitivity to λ will only depend on this particular configuration. However, a contact diagram, a t-channel diagram and a u-channel diagram have to be taken into account too as shown in Fig. 1, in which we display all the possible tree level contributions to the mentioned scattering processes in the Unitary gauge. Each of these diagrams has its own energy dependence and its own relative size, so they participate differently in the to- . This can be seen in Eqs. (1)-(4), where we show the amplitude of each diagram of the process W + W − → H H , A d , with d = s, c, t, u from s, contact, t and u channels respectively, computed consistently in the Unitary gauge: Here, g is the EW coupling constant, m W is the mass of the W boson, and s, t and u are the usual Mandelstam variables. The amplitudes for the ZZ → H H case are identical except for a global factor 1/c 2 w (with c w = cos θ w and with θ w being the weak angle), that has to be included in each amplitude, and the substitution of m 2 W by m 2 Z in the t and u channel expressions. On the other hand, the contribution of each polarization state of the initial EW gauge bosons behaves differently, not only energetically, but also in what concerns to the sensitivity to λ. There are only two polarization channels that do depend on λ: the purely longitudinal, V L V L , and the purely transverse in which both vector bosons have the same polarization, V T + V T + and V T − V T − . All the other channels have vanishing s-channel contributions and will not actively participate, therefore, in the study of the Higgs trilinear coupling, although all polarization states contribute to the total cross section. Moreover, this total cross section is dominated, specially at high energies, by the purely longitudinal V L V L configuration, and so is each diagram contribution. All these features can be seen in Fig. 2, where we display the predictions for the cross sections of W + W − → H H and ZZ → H H as a function of the center of mass energy for three different values of λ separated by polarizations of the gauge bosons, including, also, the unpolarized cross section. In this figure two things are manifest: the first one is that the V L V T configuration is indeed independent of λ. The second one is that the total cross section is clearly strongly dominated by the purely longitudinal contribution at all energies. This is a very interesting result, since it means that, if this process was measured, we would be being sensitive to the purely longitudinal configurations of the gauge bosons, and therefore to the heart of the self-interactions of the SM scalar sector.
The V L V L dominance can be understood through the inspection of the energy dependence of the longitudinal polarization vectors, ε V , at high energies. They are all proportional, for √ s m V , to a power of the energy over the mass, E V /m V . This leads to a behavior of the amplitudes, presented in Eqs. (2)-(4), for the contact, t and u channels respectively, proportional to s, and to a constant behavior with energy of the s-channel amplitude given in Eq. (1). Including the extra 1/s suppression factor to compute the cross section from the squared amplitude one obtains the energy dependence seen in Fig. 3, where we present the contribution of each diagram to the total cross sections of W + W − → H H and ZZ → H H in the SM, as well as the sum of the contact, t -channel and u-channel diagrams, (c + t + u), and the total cross section taking all diagrams into account. In this figure, we see clearly that the sum of the contact, t and u channels tends at high energy to a constant value. This happens because in the SM there is a cancellation of the linear terms in s among these three channels. In contrast, the s-channel contribution decreases as 1/s and is subleading numerically in the SM with respect to the other (c + t + u) contributions. It is only, at lower energies near the production threshold of two Higgs bosons, where the s-channel contribution is numerically comparable to the other channels and, in fact, a mild cancellation occurs between this s-channel and the rest (c + u + t). Therefore, the s-channel and in consequence λ, do not effectively participate in the constant behavior at high energies of the total cross section in the SM. At this point, it is worth recalling that these constant behaviors of the cross sections with energy are characteristic of VBS processes at high energies, precisely because of the above commented dominance of the longitudinal configurations.
When going beyond the SM by taking λ = λ SM , the previously described dependence with energy and the delicate cancellations commented above among the various contributing diagrams may change drastically. In fact, varying the size of the Higgs trilinear coupling could modify the relative importance of the contributing diagrams and, in particular, it could allow for the s-channel contribution to be very relevant or even dominate the scattering. This could happen not only at low energies close to the threshold of H H production, but also at larger energies, where the pattern of cancellations among diagrams could be strongly modified. This may lead to a different high energy behavior, and, hence, to a different experimental signature. The crucial point is that such a large deviation in λ with respect to the SM value is still experimentally possible, as the present bounds on the trilinear coupling are not yet very tight. The best bounds at present set κ = λ/λ SM ∈ [−8.2, 13.2] [54], so values of order 10 times the SM coupling are still allowed by LHC data. Then, if in the future the LHC could improve this sensitivity to lower values of λ it would be a formidable test of the presence of new physics beyond the SM. We will show next that this sensitivity can be indeed reached in the future by means of VBS. It is important to understand in more detail at this point the implications of setting λ to a different value than λ SM in the kinematical properties of the VBS processes we are studying here. For this purpose, we present in Fig. 4 the total cross section of the process W + W − → H H as a function of the center of mass energy √ s and the differential cross section with respect to the pseudorapidity η H of one of the final Higgs bosons (notice that the distribution with respect to the pseudorapidity of the other Higgs particle is the same) for different values of positive, vanishing and negative λ. 1 The results for ZZ → H H (not shown) are very similar to those of W + W − → H H . From this figure, it can be seen that, first and most evidently, the total cross section changes in magnitude and in energy dependence with respect to the SM one, as already announced. This happens especially near the H H production threshold, confirming that the sensitivity to deviations in λ with respect to the SM value is larger in this region. For the case of positive λ the total BSM cross section can be larger or lower than that in the SM, depending on the size of the deviations in λ with respect to λ SM , since in this case there is a destructive interference between the s channel contribution and the rest (c + t + u). In contrast, for the case of negative λ values, the sum of diagrams is always constructive and one obtains bigger cross sections than the SM one independently of the absolute value of the coupling. The details of these features will be extended when commenting the next figure. Regarding the angular dependence of the differential cross section, or correspondingly the distribution respect to η H also shown in Fig. 4, we see clearly that it also changes in the BSM scenarios respect to the SM one. We particularly learn from this figure that for central values of the Higgs pseudorapidity, concretely for |η H | < 2.5, it is much easier to distinguish between different values of λ. Therefore, this suggests the kind of optimal cuts in this variable η H , or the equivalent one in terms of the final particles from the Higgs decays, we should be giving to enhance the sensitivity to the signal when moving to the realistic case of the pp collisions at the LHC.
In Fig. 5 we display our predictions for the total cross section of the two relevant VBS processes as a function of κ for four different values of fixed center of mass energy √ s =  260, 500, 1000, 3000 GeV. We also display the parabolic fits that allow us to describe each of the curves to have a more analytical insight into the details of how the above commented cancellations among diagrams do actually occur. The formulas of the fits in this figure manifest that, in general, the cross section has a quadratic, a constant and a linear term in κ, coming, respectively, from the s-channel contribution, from the (c + t + u) contribution and from the interference of these two. The sign of the interference is negative for positive values of κ and positive for negative values of κ. This destructive interference for λ > 0 produces that the minima of these lines are placed at λ > λ SM . Besides, depending on the energy and on the size of κ, the behavior of the cross section will be dominantly constant, linear or quadratic in λ, and therefore the sensitivity to λ will vary accordingly. Near the production threshold, i.e., at energies around 250 GeV, two issues can be seen. The first one is that, as we already saw in Fig. 4, the differences in the cross section when we vary λ are maximal, and so will be the sensitivity to differences in this coupling. The second one is that, at these low energies, the SM, corresponding to κ = 1, suffers, as already said, a mild cancellation between the linear and the constant terms, and therefore the sensitivity to λ will be mainly quadratic. We can also see that the minima of the parabolas soften, in the sense that the variations in the cross section when we vary λ become smaller, and that their position moves from λ/λ SM close to 2 to larger values as the energy is increased. Because of this, the bigger the energy, the bigger the value of λ that maximizes the cancellations. Thus, as a first conclusion at this point, we will have to keep in mind, once we perform the full collider analysis, that the sensitivity to different values of the trilinear coupling and the issue of delicate cancellations among diagrams in VBS are clearly correlated and this will affect the final results at the LHC.
A final comment has to be made in this section, and it is that of a potential unitarity violation problem for large |λ| values in the processes of our interest here, V V → H H . To check this unitarity issue, we have evaluated the partial waves a J of the dominant polarization channels for this VBS, which, as we have said, are the longitudinal ones, i.e., V L V L → H H . These a J of fixed angular momentum J are evaluated as usual, by computing: where P J (cos θ) are the Legendre polynomials. For a given energy, √ s, we then define the unitarity violation limit as the value of λ for which |a J (s)| = 1. By doing this exercise, we find that all the partial waves |a J | that we have computed are below 0.1 for values of λ between −10 and 10 times the SM value at all energies. So, for the present study, we are safe from unitarity violation problems. For completeness, we have also made a fast estimate of the value of λ that would be required to violate unitarity in this process. For large values of |λ|, the dominant contribution to the total amplitude comes from the s-channel. This contribution, as we mentioned before, behaves, at high energies and for the purely longitudinal case, as a constant. In particular, one obtains that A s (V L V L → H H ) ∼ 6 λ for √ s m H . With this amplitude, one can compute the value of λ for which the biggest partial wave (in this case we have checked that it is the one corresponding to J = 0) becomes one. We obtain λ unit ∼ 17. Notice that this upper limit of λ is above the perturbativity limit given naively by λ pert ∼ √ 16π ∼ 7. With all these features in mind, we can move on from the subprocess level to the full process at the LHC to study the sensitivity of this collider to the Higgs self coupling in VBS processes.

Sensitivity to the Higgs self-coupling at the LHC
Once we have characterized completely the scattering V V → H H , it is time to explore the full process at the LHC to quantify how sensitive this machine could be to the Higgs trilinear coupling in VBS processes. At this point, we would like to stress again the fact that this double Higgs production channel, via the scattering of two EW gauge bosons, has been poorly studied previously in the literature, due to the fact that it provides less statistics than the GGF one. Nevertheless, now that the LHC is close to reach its nominal energy, √ s = 14 TeV, and that it is already achieving high integrated luminosities, close to L = 40 fb −1 , the possibility of measuring VBS processes, that were inaccessible before, opens up. In fact, several VBS measurements have been already performed at this collider by ATLAS [60][61][62][63][64][65][66] and CMS [67][68][69][70][71][72][73][74]. Taking this into account, and the fact that the kinematics of the VBS processes are incredibly characteristic and allow for a very efficient signal selection and background rejection, a dedicated study of the sensitivity to λ via VBS processes is on demand. This is precisely the aim of this section, in which we first promote the analysis of Section 2 to that of its LHC signal, pp → H Hjj , so that we can fully understand its behavior and properties, and then we give more quantitative and realistic results for the sensitivity to λ once the Higgs bosons have decayed. Specifically, we will focus first on the dominant Higgs decays to bottoms, leading to the process pp → bbbbjj . This process benefits from having more statistics due to the large branching ratios involved, and, because of this, it is presumably the one that will lead to the best sensitivities. We will also present results on other channels, concretely for pp → bbγ γjj , where one of the two Higgs bosons has decayed to photons, that, despite their smaller number of events, might also provide interesting results since they suffer from less severe backgrounds.
For all computations and results of the signal events we use MadGraph5@NLO [75], setting the factorization scale to Q 2 = m 2 Z and using the set of PDF's NNPDF2.3 [76]. We have found that changing the chosen value of Q 2 does not lead to relevant changes in the signal rates. Concerning the backgrounds, all of them are simulated with the same settings and PDF's as the signal, using MadGraph5@NLO as well. For the case of the multijet QCD background in the pp → bbbbjj channel, due to its complexity, we have simulated events using both Mad-Graph5 with the previous mentioned settings and PDF's, and AlpGen [77], this time choosing and selecting the set of PDF's CTEQ5L [78]. We have found agreement between the results of these two Monte Carlos, within the provided errors, in the total normalization of the cross section with the basic cuts, and in the shape of the relevant distributions. All results are presented for a center of mass energy of √ s = 14 TeV. Our study is aimed to be a first and simple approach to the sensitivity to λ in VBS processes at the LHC. This means that, in order to simplify the procedure, the analysis is done at the parton level, and no hadronization or detector response simulation are performed, leaving always room for more expert improvement towards a full and dedicated experimental study.

Study and characterization of pp → H Hjj signal events
In order to be able to estimate the sensitivity to the Higgs self-coupling in VBS at the LHC, we need to understand how the results of the previous section translate into the full process when we start with protons as initial particles. This full process, pp → H Hjj , can be produced via many different channels, and not only in VBS configurations. In fact, it is well known that this VBS subset of diagrams contributing to q 1 q 2 → q 3 q 4 H H is not gauge invariant by itself and all kinds of contributing diagrams have to be included to get gauge invariant result. This is indeed what we are doing here, since when we use MadGraph to compute the signal all kind of diagrams are included.
The crucial point regarding the phenomenological interest of VBS, that indeed motivates this work, is that the specific VBS configuration can be very efficiently selected by choosing the appropriate kinematic regions of the two extra jets variables, as it is well known [45,55,56,58]. In particular, at the LHC, the VBS topologies are characterized by large separations in pseudorapidity of the jets, | η jj | = |η j 1 − η j 2 |, and by large invariant masses of the dijet system, M jj . Imposing proper cuts over these two variables makes possible to obtain events that come dominantly from VBS processes and, as we will see later on, also to reject many background events.
The VBS processes involved in pp → H Hjj can be seen schematically in Fig. 6, where the green blob represents all diagrams in Fig. 1, including the presence of the s-channel with the generic Higgs trilinear coupling λ. This kind of processes will inherit the properties of the sub-scatterings we have studied, but will also have differences with respect to them due to the fact that we now have protons in the initial state. Then, it is important to know at this stage how close to the "pure" VBS configuration our pp → H Hjj signal is. To this end, we have generated with MadGraph5 pp → H Hjj signal events for this process for different values of λ with a set of basic cuts that allow for the detection of the final particles, given by: where p T j is the transverse momentum of the jets, η j,H is the pseudorapidity of the jets or of the Higgs bosons, and R jj is the angular separation between two jets defined as R jj = η 2 jj + φ 2 jj , with η jj and φ jj being the angular separation in the longitudinal and transverse planes, respectively. With these generated events, we have studied some relevant distributions for the signal cross section that we have found give the most efficient access to the VBS configuration in pp → H Hjj events: distributions with M jj , η jj and M H H .
In Fig. 7 we present the predictions for the cross section of the process pp → H Hjj for different values of λ as a function of the separation in pseudorapidity of the final jets | η jj | and as a function of the invariant mass of these two jets M jj . In these plots we can see that our signal is indeed dominated by the VBS configuration, since a very large fraction of the events populate the kinematic regions that correspond to VBS topologies. To have a quantitative estimation, we can take, for instance, the VBS selection cuts proposed in [58] and impose them to the events shown in Fig. 7. Thus, by imposing these cuts: we obtain that between 50% and 75% (depending on the value of λ, with closer values to 75% for the larger values of |λ|) of the events are accepted within them, which means that the VBS topologies amount, 2 at least, to half of the total cross section of pp → H Hjj . This is indeed a very interesting result, since, as we will see in the forthcoming section, the VBS cuts allow us to reduce some backgrounds even in two orders of magnitude. The fact that the signal is practically left unaffected by these cuts is an excellent outcome as the signal to background ratio will favor a better sensitivity to λ. Furthermore, knowing that the process of our interest at the LHC has a dominant VBS configuration, we would expect the translation from the subprocess results to the complete ones at this level to be straightforward. This appears to be the case, as shown in Fig. 8, where we display the predictions for the total cross section of the process pp → H Hjj as a function of the invariant mass of the diHiggs system, M H H , for different values of the Higgs self-coupling after imposing the cuts given in Eqs. (6) and (7). In these plots, it is manifest that the curves follow the same tendency as the subprocess when we vary λ. Near the H H production threshold the difference in the cross sections for different values of the coupling is more pronounced, and one can see again that the cancellations play a role in the same way we learnt at the subprocess level. The SM cross section (κ = 1, in red) lies between the κ = 0 (in green) one, which is bigger, and the Fig. 9. Predictions for the total cross section (left panel) and for the ratio of the total cross section over its SM value (right panel) as a function of the Higgs self-coupling λ with and without imposing the VBS selection cuts given in Eq. (7). Cuts in Eq. (6) have been applied and the center of mass energy has been set to √ s = 14 TeV. κ = 2 (in light blue) one, which is smaller. Again, for negative values of κ the cross section is always larger than the SM one, so we will have, for the same absolute value of the coupling, better sensitivities for negative λ values.
The issue of the cancellations that take place between the diagram that depends on λ and the rest is shown in more detail in Fig. 9. In this figure, we present the predictions for the total cross section for pp → H Hjj , and for the ratio of the total cross section over its SM value as a function of the Higgs self-coupling. We also compare the results with and without imposing the VBS cuts given in Eq. (7) to explore how the cancellation happens at the LHC, and how it depends on the selection of the VBS topologies. We learn again, that, for the same absolute value of λ, negative values give rise to larger cross sections, and therefore to better sensitivities. The smallest cross section corresponds roughly to κ ∼ 1.6, which is the value that will be harder to reach at the LHC. One may notice that this value does not coincide exactly with that in Fig. 5, even for the dominant contribution close to the threshold. This slight displacement of the minimum is due to the fact that many different topologies in addition to those of VBS contribute to this final state, in contrast with the results in Fig. 5 that took into account only VBS configurations. In fact, once we apply the VBS cuts the minimum gets closer to that of Fig. 5. Besides, and interestingly, the effect of imposing the VBS selection cuts can ameliorate the sensitivity to λ. Although the cross sections reduce in value after applying the cuts, the ratio of the total cross section for a given trilinear coupling over the SM cross section increases when we are away from the region in which the cancellations are relevant, i.e., for κ > 3 and κ < 1.
The last issue we would like to point out in this section refers to the kinematical behavior of the VBS subsystem, that is then translated to the kinematics of the final Higgs bosons. Usually, in vector boson scattering processes at the LHC, most of the energy of the initial pp state is transmitted to the radiated EW gauge bosons. This leads, as a consequence, to a very boosted system of final H H pairs, which can be profitable to select these kind of events against backgrounds. If the final Higgs particles are very boosted, their decay products, will have, in general, small angular separations. This, together with the fact that the invariant mass of the two particles that come from the Higgs decay has to lie near the Higgs mass, will allow us to characterize very efficiently the Higgs boson candidates as we will see in the next section. With this and the VBS topologies under control, we can study the full processes in which the Higgs bosons have decayed, and compute the sensitivities to λ in these realistic BSM scenarios.

Analysis after Higgs boson decays: sensitivity to λ in pp → bbbbjj events
As previously mentioned, once we have fully characterized our most basic process, pp → H Hjj , we need to take into account the Higgs decays to perform a realistic analysis at the LHC. The channel we are going to focus on is pp → bbbbjj , since the decay of the Higgs boson to a bottom-antibottom pair benefits from the biggest branching ratio, BR(H → bb) ∼ 60%. Because of this, we will obtain the largest possible rates for our signal, which will allow us to probe the broadest interval of deviations in the Higgs self-coupling.
Although this process is really interesting because of its large statistics, it is important to mention that it also suffers from having a severe background: the one coming from pure multijet QCD events. This QCD background, of O(α 3 S ) at the amplitude level, leads to the same final state as our signal, pp → bbbbjj , and, although in general they have very different kinematics, their rates are so high that some of the events can mimic the signal coming from the decay of two Higgs particles. For this reason, we need to be very efficient when applying selection cuts and criteria to be able to reject this particular background.
We learnt in the previous sections that our signal is very dominated by the VBS configuration. Oppositely, the multijet QCD background is composed primarily by topologies that do not share kinematical properties with VBS processes. This is the reason why we will first select those QCD events that can be misidentified as those signal events coming from VBS, and take them as a starting point to perform our more refined study of the signal and background.
To have a first insight on how efficient the VBS selection criteria are, we have generated with MadGraph5 ten thousand events for our signal, pp → H Hjj → bbbbjj in the SM, i.e., κ = 1, and for the multijet QCD background with a set of basic cuts that ensure the detection of the final state particles: where p T j ,b is the transverse momentum of the jets and bottoms, η j,b are the pseudorapidities of the jets or of the bottom particles, and R ij is the angular separation between the i and j particles.
In Fig. 10 we display the localization of these events in the | η jj | − M jj plane, the two variables that better characterize the VBS processes. One can see, indeed, that the QCD events populate mostly the region of small invariant masses of the dijet system and of small differences in pseudorapidity of the jets, as opposed, precisely, to the signal events. Thus, imposing the proper VBS cuts, like those in Eq. (7), should relevantly reduce the QCD background leaving the signal nearly unaffected.
In Fig. 11 we aim precisely to see this effect, since we present the same set of events as in Fig. 10 for the QCD background and for the signal highlighting in orange those events that fulfill the VBS selection criteria given in Eq. (7) as an example. This time we show the results in the M bb 1 − M bb 2 plane, where M bb 1,2 are the corresponding invariant masses of the two bottom pairs that are the best candidates to come from the decay of a Higgs boson, as we will see later.
The first thing one can observe in both plots of Fig. 11 is that very few QCD events survive the imposition of the VBS cuts, whereas practically all events of the signal do. The concrete fraction of the events (A) that survive in both cases is also presented in the plots. We call A VBS the acceptance of the VBS cuts, defined as i.e., the ratio between the cross section of the process after applying the VBS cuts like those in Eq. (7) over the cross section of the process without having applied them. The basic cuts are imposed in all cases. Taking a look at these numbers, we see that 60% of the signal events pass these cuts while only 9% of the QCD events do. At this point, one might wonder whether these results are very dependent on the specific VBS cuts we impose or not. In Table 1 we show the predictions for the acceptances, A VBS , of different sets of VBS selection cuts, i.e., different cuts in | η jj | and in M jj , for both the multijet QCD background and the signal with κ = 1. From those predictions we can see that all the sets of cuts considered lead to very similar results: around 60% of the signal fulfills the VBS selection criteria whereas a 5-10% of the multijet QCD background does. We have checked that for other values of κ the acceptance for the signal varies between a 55% and a 75%. From now on we will apply the VBS selection cuts given in Eq. (7), since this set is well explored in the literature and qualitatively provides the same results as the other sets of cuts that we have analyzed.
The second issue that we can notice about Fig. 11 is that, again, the QCD events populate a very different region of this plane than those of the signal. QCD events tend to lie on low values of M bb i , somehow away from the region close to the [M bb 1 = m H , M bb 2 = m H ] point in the M bb 1 − M bb 2 plane, in which most of our signal settles. Evidently, two particles coming from the decay of a Higgs boson should have a total invariant mass value near the Higgs boson mass, m H , as our signal does. This motivates the next selection criteria we are going to apply, following the search strategies of ATLAS [52] and CMS [50] for double Higgs production, that are aimed to efficiently identify the H H candidates.
The H H candidate identification criteria are also based on what we have learned in the previous sections. Logically, each H candidate corresponds to a b-quark pair, and therefore we first need to define how we are going to pair the final b-quarks. From now on, it is worth mentioning that we will not distinguish between bottom and anti-bottom, similarly to what is done in experimental analyses. Therefore, with four bottom-like particles in the final state we have three possible double pairings. Among these three possibilities, we select the one in which the values Orange dots correspond to those events that pass the implemented VBS selection cuts given in Eq. (7). Cuts in Eq. (8) have been implemented. The value of the acceptance A is also included. The red cross represents the value of the Higgs mass The center of mass energy has been set to √ s =14 TeV.
where the super-indices l and s denote, respectively, leading and subleading, defining the leading b-quark pair as the one with largest scalar sum of p T . One might notice that the requirement of small angular separation between the two b-quarks of a pair, and the fact that the invariant mass of each b-quark pair has to lie near the mass of the Higgs, are encoded in the ˆ R bb and in the χ H H cuts, respectively. The latter is equivalent to impose that the events in the M bb 1  Nevertheless, although multijet QCD events represent the most severe background, there are other processes that can fake our signal. One of them is the tt background, with the subsequent decays of the top quarks and W bosons, tt → bW +b W − → bbbbjj . This is, however, a very controlled background, since it is well suppressed by non-diagonal CKM matrix elements and its kinematics are radically different than those of VBS. Starting from a cross section of 5.4 · 10 −5 pb with all the basic cuts in Eq. (8) applied, one ends up in 1.7 · 10 −7 pb after applying the H H candidate cuts, and in 2.0 · 10 −10 pb after applying the VBS cuts afterwards. Therefore, since this background is five orders of magnitude smaller than the smallest of our signals, we will neglect it from now onwards. Finally, we still have to deal with other potentially important backgrounds corresponding to pp → H Zjj → bbbbjj and pp → ZZjj → bbbbjj . These two H Z and ZZ production processes, receiving contributions of order (α · α S ) and (α 2 ) at the amplitude level, also drive to the same final state as our signal and may give rise to similar kinematics, since they can also take place through VBS configurations. In fact, their rates are very close to those of our signal after applying the VBS selection cuts, that reduce these backgrounds less efficiently than the multijet QCD one. However, we can again take advantage of the fact that the b-quark pairs have to come from a Higgs boson with a well defined mass. Therefore the H H candidate cuts should allow us to reject these backgrounds.
In Table 2 we present the cross sections of the multijet QCD background, of the combined pp → H Zjj → bbbbjj and pp → ZZjj → bbbbjj background and of the signal with κ = 1, with the basic cuts already set, after applying each of the cuts in Eqs. (10)- (13) subsequently. This way, we see the reduction factor after each cut, and the total cross section of both signal and background once we have performed our complete H H candidate selection. We show as well the effect of applying the VBS cuts given in Eq. (7) afterwards, since we have checked that both sets of cuts (H H candidate cuts and VBS cuts) are practically independent. Thus, we have the total cross sections of the two main backgrounds and of our SM signal after applying all the selection criteria. In Table 3 we provide the total cross sections of the signal for all the values of λ considered in this work, again after applying all the selection criteria, for comparison.
From the results in Table 2 we can learn that the sum of the two backgrounds, ZHjj +ZZjj , is under control after applying the H H candidate cuts, since its cross section lies an order of Table 2 Predictions for the total cross section of the multijet QCD background, of the combined pp → H Zjj → bbbbjj and pp → ZZjj → bbbbjj background and of the signal with κ = 1 after imposing each of the cuts given in Eq. (8) and in Eqs. (10)-(13) subsequently. We show as well the total cross section after applying, afterwards, the VBS selection cuts in Eq. (7).  (13) 7 .9 · 10 −2 8.6 · 10 −6 9.0 · 10 −5 VBS cuts in Eq. (7) 6 .8 · 10 −3 5.5 · 10 −6 4.1 · 10 −5 Table 3 Predictions for the total cross section of the signal pp → bbbbjj after imposing all the selection criteria, VBS cuts given in Eq. (7) and H H candidate cuts given in Eqs.  Table 4 Predictions for acceptances of the VBS cuts given in Eq. (7) magnitude below the SM signal. On the other hand, the multijet QCD background remains being very relevant even after imposing all the selection criteria. However, as we will see later, the total reduction that it suffers still allows to be sensitive to interesting values of κ even for low luminosities. This reduction, along with that suffered by the ZHjj +ZZjj backgrounds and with that suffered by the SM signal, is presented in Table 4. There we show the acceptances of the VBS cuts and the H H candidate cuts separately and together for the multijet QCD background, for the combined pp → H Zjj → bbbbjj and pp → ZZjj → bbbbjj background and for the SM signal, for comparison. It must be noticed that other backgrounds apart from those having the same final particle content as our signal can contribute relevantly. This would be the case if some final state particles were misidentified, leading to a "fake" bbbbjj state. The most important of these backgrounds is the production of a tt pair decaying into two b quarks and four light jets, tt → bW +b W − → bbjjjj with two of these final light jets being misidentified as two b jets. In order to estimate the contribution of this background, we have generated with MadGraph5 tt → bW +b W − → bbjjjj events applying first the minimal cuts |p T j,b | > 20 GeV, |η j,b | < 5 and R jj,bj,bb > 0.2, with a total cross section of 246 pb. Applying a mistagging rate of 1% per each light jet misidentified as a b jet, we obtain 246 · (0.01) 2 = 2.5 · 10 −2 pb as starting point to compare to our main multijet QCD bbbbjj background. Now we need to apply our selection cuts described in Eqs. (7)- (13) to see their impact on this particular background. We apply first the VBS selection cuts asking that there is at least one pair of light jets fulfilling the criteria in Eq. (7). These cuts reduce the cross section to 1.3 · 10 −4 pb. Now analyzing the events that pass the VBS cuts, if there is only one pair of "VBS-like" light jets, the other two light jets are identified as b quarks. If there is more than one, we select as b quarks those that minimize |M pp 1 − M pp 2 |, with p = b, j among all possibilities. Once we have characterized our two light jet candidates and our four b-quark candidates, we proceed with the H H candidate selection cuts. This way, applying subsequently the criteria explained in Eqs. (10)-(13), we obtain the following cross sections: 1.2 · 10 −5 pb (p T b ), 2.5 · 10 −6 pb (ˆ R bb ), 4.4 · 10 −7 pb (p T bb ) and finally 2.1 · 10 −8 pb (χ H H ). Therefore, since this tt background is five orders of magnitude below our main considered background, whose final cross section given in Table 2 is 6.8 · 10 −3 pb, we conclude that it can be safely neglected.
We have also considered the possible backgrounds coming from multijet QCD processes leading to different final states than that of bbbbjj , such as 6j and bbjjjj , in which some of the final state light jets are again misidentified as b jets. To estimate their contribution to the background we have used the total cross sections of these processes given in [77]. These are, for a center of mass of 14 TeV, 1.3 · 10 5 pb and 7.5 · 10 3 pb, respectively. If we apply now the corresponding misidentification rates we end up with 1.3 · 10 5 · (0.01) 4 = 1.3 · 10 −3 pb for the case in which we have six light jets, and 7.5 · 10 3 · (0.01) 2 = 7.5 · 10 −1 pb for the case in which we have two b jets and four light jets. We now assume that the selection cuts we specify in Eqs. (7)-(13) will have a similar impact on these backgrounds as they do on the multijet QCD production of four b jets and two light jets, since they all take place through similar QCD configurations. Thus, applying the corresponding acceptance factor of these cuts we obtain the following total cross sections: 1.3 · 10 −3 · 1.1 · 10 −5 = 1.4 · 10 −8 pb for the six light jets case and 7.5 · 10 −1 · 1.1 · 10 −5 = 8.2 · 10 −6 pb for the two b jets and four light jets case. Both of these cross sections are more than three orders of magnitude below that of the bbbbjj background, so we conclude that they can also be safely neglected without introducing big uncertainties.
Once we have the possible backgrounds under control, we can move on to fully explore the sensitivity to the Higgs self-coupling λ in pp → bbbbjj events. In Fig. 12 we display the predictions for the total cross section of the total SM background (the sum of multijet QCD background and the combined ZHjj + ZZjj backgrounds) and of the signal for different values of λ as a function of the invariant mass of the four-bottom system M bbbb . These distributions should be the analogous to those in Fig. 8 after the Higgs boson decays, as it is manifest since the signal curves follow the same tendency and are very similar except for the global factor of the Higgsto-bottoms branching ratio. In this figure we can also see that the total SM background is of the same order of magnitude than the κ = 10 and κ = −5 signals, and it is even below the κ = −10 signal prediction. This is a very interesting result, since it means that if, for example, the true value of λ was minus five times that of the SM, the LHC should be able to measure twice as many events as those expected from the SM background only in this VBS configuration. Similar conclusions can be extracted for other values of κ.
Given the encouraging previous results, our last step is to give quantitative predictions for the sensitivity to λ in pp → bbbbjj processes at the LHC. To this end, we compute the statistical significance S stat , as defined in [79] by: where N S and N B are the number of events of signal and background, respectively. Notice that for N S /N B 1, this definition of S stat tends to the usual N S / √ N B expression. This computation is going to be performed for four different values of the luminosity: L = 50, 300, 1000, 3000 fb −1 , that correspond to a near-future LHC value for the current run (50 fb −1 ), and to planned luminosities for the third run (300 fb −1 ) and the High-Luminosity LHC (HL-LHC) (1000 and 3000 fb −1 ) [80].
In Fig. 13 we present the results of the statistical significance of our signal, S stat , in pp → bbbbjj events as a function of the value of κ, for the four luminosities considered. We display as well a closer look for the values of κ ranging between 0.5 and 2.5, interesting for an elevated number of well motivated BSM models. In the lower part of the left panel we also present the corresponding predictions for the total number of signal events, N S , as a function of κ. The marked points correspond to our evaluated predictions. We show as well, in the right panel of this figure, our predictions for the value of the total integrated luminosity, L, as a function of the value of κ as well, that will be required to obtain a sensitivity to a given κ in pp → bbbbjj events at the 3σ and 5σ level. In this plot, we have also marked the areas in luminosity where the number of predicted signal events N S is below 1, 10 and 100, respectively, to get a reference of the statistics obtained.
From these plots, we can extract directly the conclusions on the sensitivity to λ in VBS processes at the LHC in pp → bbbbjj events. The first thing one might observe is the high statistics and significances of the signal for most of the studied cases, except for the region close to the SM value, say for κ between 1 and 2. Studying carefully this particular region of the parameter space, we conclude that it is the most challenging one to access at the LHC, since all the predicted statistical significances given for κ ∈ [0.5, 2] are below 2σ even for the highest luminosity considered. The second one is that, for the same absolute value of the coupling, the sensitivities to negative values of κ are higher than to positive values of κ. The third conclusion is that the LHC should be sensitive to very broad intervals of κ, even for the lowest luminosity considered, L = 50 fb −1 , with high statistical significance. These means that VBS processes could allow us to probe the value of λ with very good accuracy in the near future. More specifically, in Table 5 we show the summary of the predictions for the values of κ ≡ λ/λ SM that the LHC would be able to probe in pp → bbbbjj events, with a sensitivity equal or better than 3σ (5σ ) for the four luminosities considered: L = 50, 300, 1000, 3000 fb −1 . These results are indeed very interesting, since the sensitivities to λ that one can obtain from studying VBS double Higgs production are very promising even for the lowest luminosity considered 50 fb −1 . The ranges of λ that the LHC could be able to probe in this kind of processes indicate that it is worth to study VBS as a viable and useful production mechanism to measure the Higgs trilinear coupling. On the other hand, it can be seen that the HL-LHC should be able to test very small deviations in the value of the Higgs self-coupling and that it should be sensitive to all the explored negative values for κ. Although the present work is a naive study, since it is performed at the parton level and does not take into account hadronization and detector response simulation, the results in Table 5 show that the VBS production channel could be very promising to measure the true value of λ, and, therefore, to understand the nature of the Higgs mechanism.

Analysis after Higgs boson decays: sensitivity to λ in pp → bbγ γjj events
The pp → bbbbjj process is, as we have seen, a very promising channel to study the Higgs self-coupling at the LHC due to its large event rates. However, it is clear that it suffers from quite severe backgrounds, coming specially from multijet QCD events, so one could think of studying complementary channels with smaller rates but with a cleaner experimental signature. This is the reason why we would like to explore the case in which one of the Higgs bosons decays to a b-quark pair, as before, while the other one decays to two photons through gauge bosons Table 6 Predictions for the total cross section of the signal pp → bbγ γjj after imposing all the selection criteria, VBS cuts given in Eq. (7) and cuts given in Eqs. (16) and (17) for all the values of κ considered in this work: κ = 0, ±1, ±2, ±5, ±10. The cross section of the SM background for this same cuts amounts to σ Background = 1.4 · 10 −6 pb. Basic cuts in Eq. (15) are also applied. and fermion loops. This implies a large reduction factor in statistics due to the comparative low branching ratio BR(H → γ γ ) ∼ 0.2%, a factor 0.003 smaller than that of H → bb. The analysis of the process pp → bbγ γjj implies to go through its main backgrounds as well. We will consider in this section the same background ZH of the previous case, since the ZH final state can also lead to processes with two photons and two bottoms, pp → H Zjj → bbγ γjj , coming from the decays of the H and the Z. In addition, we also consider the mixed QCD-EW pp → bbγ γjj background, of O(α · α 2 S ) at the amplitude level, that should be the most severe one.
As we did before, to study signal and background, we first need to establish a set of cuts that ensure particle detection, so we apply the following basic cuts: and afterwards, to reject the QCD-EW and the ZHjj backgrounds we will apply first the VBS cuts given in Eq. (7) and subsequently the following kinematical cuts given by CMS in [53]: where l and s stand for leading (highest p T value) and subleading photons, and where M γ γ is the invariant mass of the photon pair. The final ingredient is to apply the χ H H cut, taking now into account that we have a b-quark pair and a photon pair in the final state: This ensures that the two b-quarks and the two photons come from the decay of a Higgs particle. Once again, there might be important background contributions from multijet QCD processes leading to different final states than that of bbγ γjj , such as 6j and bbjjjj , in which some of the final state light jets are again misidentified as b jets and some are misidentified as photons. Taking again as the presumably leading QCD background processes the production of six light jets and of two b jets and four light jets, applying a misidentification rate of 0.1% per each jet misidentified as a photon, and considering a similar reduction factor after our selection cuts as before, since the selection cuts are very similar, we obtain: 1.3 · 10 5 · (0.01) 2 · (0.001) 2 · 1.1 · 10 −5 = 1.4 · 10 −10 pb for the six light jets case and 7.5 · 10 3 · (0.001) 2 · 1.1 · 10 −5 = 8.2 · 10 −8 pb for the 2b4j case. Again in both cases the final cross sections are more than one order of magnitude smaller than the main background we have considered, being of order 10 −6 pb, concluding again that they can be neglected as well.
Having all this in mind, we present in Fig. 14 the predictions for the total cross section of the process pp → bbγ γjj as a function of the invariant mass of the bbγ γ system M bbγ γ , for different values of the Higgs self-coupling λ. We also display the prediction for the total SM  Table 7 Predictions for the values of κ ≡ λ/λ SM that the LHC would be able to probe in pp → bbγ γjj events, with a sensitivity equal or better than 3σ (5σ ) for the four luminosities considered: L = 50, 300, 1000, 3000 fb −1 .
background (sum of the QCD-EW and the ZHjj background) for comparison. Once again, one can see that the signal distributions for different values of κ are very similar to those shown in Fig. 8, and that the main difference is due to the reduction factor of the branching ratios into photons and into b-quarks. They are very similar, too, to the results of the bbbbjj final state, in Fig. 12, although two-three orders of magnitude smaller. The background is, however, very different with respect to the one for bbbbjj events. It is smaller in comparison with the signal, specially at high M bbγ γ , since it decreases much more steeply. Therefore, we would expect to have good sensitivities to the Higgs self-coupling despite the lower rates of the process involving photons. For completeness, we display in Table 6 the predictions for the total cross section of the signal, for the set of κ values considered, and after applying all cuts given in Eq. (7) and in Eqs. (15)- (17). The prediction for the cross section of the total SM background for this same cuts amounts to σ Background = 1.4 · 10 −6 pb. In Fig. 15 we show the predictions for the statistical significance S stat , computed in the same way as in the previous section, making use of Eq. (14), for the four luminosities considered previously, L = 50, 300, 1000, 3000 fb −1 and taking again a closer look for the values of κ ranging between 0.5 and 2.5. We also show the predictions of the final number of signal events, N S as a function of κ, for these same luminosities. On the right panel of this figure we present the prediction for the value of the luminosity that will be required to probe a given κ value with sensitivities at 3σ and 5σ , as a function of the value of κ. In these plots, due to the lower statistics of this process, some of the computed significances correspond to scenarios in which there is not even one signal event. The concrete predictions for these signal event rates can be read from the lower plot of the left panel. Fig. 15. Prediction of the statistical significance, S stat , of the process pp → bbγ γjj for the four luminosities considered L = 50, 300, 1000, 3000 fb −1 (left panel) and of the value of the luminosity that will be required to probe a given κ at the LHC at 3σ and at 5σ (right panel), as a function of the value of κ. The marked points represent our evaluations. In the left panel, a zoom is performed on the interesting values of κ ranging between 0.5 and 2.5. The shadowed areas in the right panel correspond to the regions where the number of predicted signal events N S is below 1, and 10. The center of mass energy has been set to √ s = 14 TeV.
Taking a look at these figures, we can again extract the conclusions on the sensitivity to the Higgs self-coupling at the LHC in VBS processes, this time in pp → bbγ γjj events. One might notice that, although the results are less encouraging than those of pp → bbbbjj events, this channel could also be very useful to measure the value of λ. Analogously to the previous section, in Table 7 we present the values of κ ≡ λ/λ SM that would be accessible at the LHC in these type of events, pp → bbγ γjj , with a statistical significance equal or better than 3σ (5σ ), for the four luminosities considered.
These results show again that the values of κ that can be probed in the future at LHC through the study of VBS processes leading to the final state bbγ γjj could be very competitive as well. Except for the lowest luminosity considered, L = 50 fb −1 , where the signal rates found at the parton level are too low as to survive the extra factors suppression due to the missing detector efficiencies, hadronization effects etc, the sensitivities found point towards the potential of VBS processes in order to obtain a precise measurement of λ. The values close to the SM value, are, again, very challenging to reach at the LHC, since the statistical significances of κ ∈ [0.5, 2] are always below 2σ for this case as well. However, the HL-LHC should be able to probe deviations in λ very efficiently in this channel.

Discussion
Finally, to close this section of results, we find pertinent to discuss on how the precision of our predictions could be improved by including additional considerations. We comment here just on those that we consider are the most relevant ones.
1.-Our computation of the H Hjj signal rates incorporates just those coming from the subprocess qq → H Hjj , which includes VBS, but this is not the only contributing channel. It is well known that also the subprocess gg → H Hjj , initiated by gluons, does contribute to these signal rates, and it is also sensitive to large BSM λ values [24]. Although it is a one-loop subprocess, mediated mainly by top quark loops, it provides a sizable contribution to the total H Hjj signal cross section. For instance, for the case of λ = λ SM , the total cross section at the LHC with √ s = 14 TeV is, according to [33], 5.5 fb from gg → H Hjj to be compared with 2 fb from VBS. Therefore, when considering both contributions to the signal, the sensitivity to λ presumably increases. However, we have explicitly checked that once we apply our optimized VBS selection cuts summarized in Eq. (7) and in Table 1, we get a notably reduced cross section for this gg subprocess. In particular, our estimate of the signal rates at the LHC with √ s = 14 TeV from gg → H Hjj , after applying the stringent M jj > 500 GeV cut and using the results in [33] for the M jj distribution, gives a strong reduction in the corresponding cross section, and leads to smaller rates for gg than those from VBS by about a factor of 20. Therefore its contribution to the signal rates studied here can be safely neglected, and no much better precision will be obtained by including this new contribution in the signal rates. We have also checked that this finding is true for other BSM values of λ.
2.-When considering next to leading order (NLO) QCD corrections in our estimates of both the signal and background rates, we expect some modifications in our results. These can be very easily estimated, as usual, by using the corresponding K-factors. Thus, for instance, for the leading bbbbjj final state, we can include these NLO corrections by taking into account the K-factors for the VBS signal and for the main background from multijet QCD. For the signal we take the K-factor from [26], given by K VBS = 1.09. For the QCD-multijet background the corresponding K-factor is, to our knowledge, not available in the literature, and different choices are usually assumed. We consider here two choices: K QCD = 1.5, and another more conservative one K QCD = 3. This implies that our predictions for the signal rates are practically unchanged, but those for the background rates are increased by a factor of 1.5 and 3 respectively. This modifies our predictions for the statistical significance of the bbbbjj signal, from the S stat results given Fig. 13 to S NLO stat ∼ K VBS / K QCD S stat ∼ 0.9 S stat and 0.6 S stat for K QCD = 1.5 and K QCD = 3, respectively. For instance, for the high luminosity considered of 1000 fb −1 we get sensitivities of κ > 3.8(4.3) for K QCD = 1.5, and of κ > 4.5(4.8) for K QCD = 3, both at the 3σ (5σ ) level, to be compared with our benchmark result in Table 5 of κ > 3.7(4.2). Therefore ignoring these NLO corrections does not provide large uncertainties either. 3.-When including b-tagging efficiencies in our estimates of the bbbbjj signal and background rates, our predictions of the statistical significance do also change. However, an estimate of this change can be easily done by adding the corresponding modifying factors. For instance, by assuming well known b-tagging efficiencies of 70%, that apply to both the signal and background, the two rates are reduced by a factor of 0.7 4 ∼ 0.24. Therefore we get a reduced statistical significance of S b−tag stat ∼ 0.24/ √ 0.24 S stat ∼ 0.5 S stat with respect to the ones that we have reported previously. This factor of 0.5 will change our predicted sensitivities to BSM λ values. Again, as an example, for the considered luminosity of 1000 fb −1 we get sensitivities of κ > 4.3(4.9) at the 3σ (5σ ) level to be compared with our benchmark result in Table 5 of κ > 3.7 (4.2). Similarly, considering also photon-identification efficiencies (also called in this work γ -tagging) of 95%, as presented in the literature, we get reduced signal and background rates for the bbγ γjj final state by a factor of 0.7 2 × 0.95 2 ∼ 0.44. Accordingly, we obtain a reduction in the statistical significance of the bbγ γjj events, given by S b,γ −tag stat ∼ 0.44/ √ 0.44 S stat ∼ 0.7 S stat with respect to our results reported in the pages above. The changes in the sensitivities to κ can be easily derived. using the same illustrative example, for 1000 fb −1 of luminosity, we get sensitivities of κ > 6.0(8.0) at the 3σ (5σ ) level to be compared with our benchmark result in Table 7 of κ > 4.6(6.0). 4 with E being the energy resolution, which in this case was set to 0.05, leading to a mass resolution of 5% of the Higgs mass, since this value optimizes the selection efficiency and could be useful for future experiments with better energy resolution. However a more realistic choice, given the current energy resolution at the LHC experiments, could rather be E · m H = 0.1 × m H GeV ∼ 12.5 GeV. We have redone the analysis with this alternative and more conservative choice and we have obtained, as expected, a reduced statistical significance. The signal rates do not change (we still get 4.1 × 10 −5 pb), but the main QCD-background does (we get 1.8 × 10 −2 pb instead of our benchmark value of 6.8 × 10 −3 pb). This translates in a reduction of the significance given by a factor S The implication of this reduction can directly be seen as a modification of the sensitivity to κ. Once again, for our benchmark case of 1000 fb −1 , we obtain sensitivities of κ > 4.0(4.5) at the 3σ (5σ ) level to be compared with our benchmark result in Table 5 of κ > 3.7(4.2). Apart from redoing the analysis for this 10% resolution, 4 we have also studied other possible and realistic values such as E = 20% and E = 30%, to have a better idea of the implications of the value of the mass determination uncertainty in our predictions. The results for both of our signals are shown in Fig. 16 by the green lines and green shaded areas, where we present the values for the statistical significance at 1000 fb −1 as a function of the value of κ for different energy resolutions of E = 5% (original scenario throughout the work), 10%, 20% and 30% (the purple lines and purple areas of this figure will be discussed in the next point of this discussion section). One can see that, as expected, the statistical significance decreases as the energy resolution worsens, but in any case, from the most optimistic case ( E = 5%) to the less optimistic one ( E = 30%), we only obtain a reduction factor of at most 0.4 in the statistical significance. 5.-Another important point that might change significantly our predictions is that introduced by the Higgs mass reconstruction uncertainty coming from detector effects. To estimate this uncertainty, we have applied a Gaussian smearing to the energy of all final state partons. Following [81], this Gaussian dispersion has been introduced as 1/ √ 2πσ · e −x 2 /(2σ 2 ) , with σ = 0.05 · E j,b for the energy dispersion of the final light and b jets and with σ = 0.02 · E γ for the energy dispersion of the final photons. We have performed this for each studied signal and for their corresponding backgrounds in order to characterize the impact that these detector effects have regarding the distribution of our events on the relevant kinematical   variables. In Fig. 17  The orange points correspond to those events that fulfill the VBS selection criteria. The impact of these VBS cuts does not change appreciably after the smearing, not on the signal nor on the background events. The selection of the Higgs candidates in the case of the bbbbjj signal is performed as explained in the text, following the minimization of |M bb 1 − M bb 2 |. This is the reason why we obtain several points distributed in the diagonal of the left panel. As expected, the detector effects translate into a dispersion of the signal points from the Higgs mass point outwards. In the bbbbjj case, the dispersion is isotropic, since the smearing affects all four b-quarks in the same way, whereas in the bbbγ γjj case, the dispersion in the M bb direction is bigger with respect to that in the M γ γ direction, accordingly to the difference in the energy resolution of b-quarks and photons in the detector. These results are compatible to those obtained in reference [38]. In any case, both signals seem to lie inside a circle of radius around 12 GeV, which corresponds to a 10% of the Higgs mass value. This suggests that the effects of the smearing on our predictions of the statistical significance will severely depend on the E we use in the χ H H selection cut, and, in principle, we expect that for E = 10% we will obtain the best sensitivities. This is so because, for this E = 10%, we select the minimum possible number of background events compatible with selecting all of our signal events simultaneously.
In order to better understand the impact of the E value once the detector effects have been taken into account, we present in the purple lines and purple shaded areas of Fig. 16 the values for the statistical significance at 1000 fb −1 as a function of the value of κ for different energy resolutions of 5% (original scenario throughout the work), 10%, 20% and 30% after the smearing on the energy of all final state partons has been applied. Is it clear from this figure that, indeed, taking E = 10% in the bbbbjj case maximizes the statistical significance once the detector effects are included. In the bbγ γjj case (notice that the purple area overlaps with the green one) the E = 5% is still the value that gives the best sensitivities, since the signal to background ratio is larger. In any case, from the upper green line to the upper purple line, there is at most a reduction factor of 0.4 in the statistical significance. 6.-Considering in addition the effects from showering and clustering of the final jets will presumably change our naive parton level predictions. As we have said, it is not the purpose of this paper to provide a full complete analysis including these important effects. It is clearly beyond the scope of this work and they will require a more sophisticated and devoted analysis with full computing power and the use of additional techniques like Boost Decision Trees (BDT) and others. This is particularly involved if we wish to control efficiently the background form QCD-multijets and, consequently, we have postponed this full analysis for a future work in collaboration with our experimental colleagues. 5 However, to get a first indication of the importance of these effects in the signal rates, we have performed a computation of the bbbbjj signal MadGraph events after showering with PYTHIA8 [82] and clustering with MadAnalysis5 [83][84][85][86] by using the anti-k t algorithm with R = 0.4. We have performed this signal estimate for one BSM example with κ = 5 and we have obtained that the cross-section after applying our basic and VBS cuts is 3.0 · 10 −3 pb if we include show-ering+clustering, which should be compared with our parton level estimate of 3.7 · 10 −3 pb, i.e. without applying showering+clustering. Therefore, the effect from showering+clustering at this signal level is not very relevant. However, it is expected that it could be relevant in the H H selection candidates and, as we have said, in the reduction efficiency of the QCDmultijet background. Nevertheless this is left for our future project.
Finally, to conclude this discussion section and in order to give a more accurate and realistic prediction, all the above mentioned considerations must be taken into account simultaneously. To this aim, we present, in Fig. 18, the predictions of the statistical significance as a function of the value of κ at 1000 fb −1 for two comparative scenarios: the original analysis from LO parton level predictions (dark grey line) and the analysis performed after taking into account the main distorting effects which are the tagging efficiencies of the final state particles, described in point 3.-of this discussion, the NLO corrections described in point 2.-and the estimation of the detector effects, introduced in point 5.-with a 10% Higgs mass determination uncertainty. We give these predictions for both of the studied signals: pp → bbbbjj (left panel) and pp → bbγ γjj (right panel). The main conclusion is that the biggest uncertainty in our predictions comes from the fact that we are not taking into account, a priori, detector effects. We have already seen that this can reduce the statistical significance by a factor of 0.4. The second biggest source of uncertainty is the choice of the value of the Higgs mass resolution, E . Taking a 10% mass resolution instead of a 5% can account for a reduction of 0.7 in the statistical significance. Similarly, the b-tagging efficiencies in the bbbbjj case can lead to a similar reduction factor of 0.7. Finally, the NLO corrections play the less relevant role when estimating the uncertainties of the calculation. All the main effects together lead to a reduction factor of at most 0.2 in the statistical significance for pp → bbbbjj and of at most 0.5 for pp → bbγ γjj . The corresponding changes in the sensitivities to κ can be easily derived from Fig. 18. Using the same illustrative example, for 1000 fb −1 of luminosity, we get sensitivities to κ > 6.2(7.7) at the 3σ (5σ ) level to be compared with our benchmark result in Table 5 of κ > 3.7(4.2) for the bbbbjj case and of κ > 7.7(9.4) at the 3σ (5σ ) level to be compared with our benchmark result in Table 7 of κ > 4.6(6.0) for the bbγ γjj one.
Based on the discussion above we believe that a more dedicated analysis, including more accurately all the considerations above with showering, clustering, and detector effects, and op-timizing the selection criteria accordingly, 6 might lead to a sensitivity to the Higgs self-coupling of the same order of magnitude, although a bit smaller, than the one obtained with our naive original analysis. We believe that our findings indicate that double Higgs production via vector boson scattering is a viable and promising observable to measure the Higgs self-coupling in BSM scenarios.

Conclusions
Being able to determine with precision the value of the Higgs self-coupling λ would allow us to understand the true nature of the Higgs mechanism and, therefore, of the scalar sector of the SM. In particular, an independent measurement of λ and m H will be crucial in this understanding. At the LHC, the most sensitive channel to this coupling λ is that of double Higgs production, that can take place through several initial configurations. Most of the theoretical and experimental studies of H H production focus on gluon-gluon fusion since it benefits from the largest rates. Nevertheless, double Higgs production by vector boson scattering has important advantages with respect to gluon-gluon fusion that, despite its lower statistics, make of it a very promising and competitive channel to probe the Higgs self coupling at the LHC. These features have motivated our study here.
In the present work, we have analyzed the sensitivity to λ in double Higgs production via vector boson scattering at the LHC, taking advantage of the fact that these processes have very characteristic kinematics that allow us to select them very efficiently against competing SM backgrounds. We have first explored and characterized the VBS subprocesses of our interest, W W → H H and ZZ → H H , both for the SM with λ = λ SM , and for BSM scenarios with λ = κ λ SM , considering values of κ between −10 and 10, to move afterwards to the LHC scenario. We have then studied the process pp → H Hjj , in order to understand the properties of this scattering, and finally we have explored and provided quantitative results for the sensitivity to the Higgs self-coupling after the Higgs decays.
We have focused mainly on the pp → bbbbjj process since it benefits from the largest rates. After applying all our selection criteria, based on the VBS characteristic kinematical configuration and in the H H candidates reconstruction, we give predictions for the sensitivity to λ in pp → bbbbjj events at the parton level for √ s = 14 TeV and for different future expected luminosities: L = 50, 300, 1000, 3000 fb −1 . Our main results for this channel are summarized in Table 5 and in Fig. 13, in which we present the values of κ that the LHC would be sensitive to at the 3σ and at the 5σ level. The sensitivities we obtain here, at the parton level, even for the lowest luminosity, are very encouraging and clearly invite to explore this channel further with the new technology applied to control the QCD-multijet background, including hadronization and detector effects, which will allow us to get a fully realistic result. Furthermore, our predictions show that the HL-LHC should be able to probe small deviations in λ respect to the SM value, reaching very good sensitivities with the highest luminosity to up to κ > 3.2 (3.7) at the 3σ (5σ ) level in the best scenario for positive values. In the case of negative values, the HL-LHC would be sensitive to all the κ < 0 values that have been considered in this work.
We give as well predictions for pp → bbγ γjj events, also at the parton level and for √ s = 14 TeV, due to the fact that it provides a cleaner, although with smaller rates, signature. The results of the sensitivities to the Higgs self-coupling in this channel, after applying the proper selection criteria, are collected in Table 7 and in Fig. 15. Again, we obtain very promising results for bbγ γ events, except for the lowest luminosity considered, where the signal rates found are too low. Interestingly, we show that the statistical significance grows faster with luminosity in this channel. This would imply that for bigger luminosities, very small deviations in λ could also be measured in pp → bbγ γjj events. In particular, for the highest luminosity considered, our predictions show that the HL-LHC could reach sensitivities to κ > 3.8 (4.7) at the 3σ (5σ ) level in this channel.
Furthermore, we give predictions for the interesting case of L = 1000 fb −1 of how the sensitivity to λ will change from our naive parton level results when taking into account the main distorting effects. We have discussed the impact of b-tagging and of γ identification efficiencies, of detector effects and of Higgs mass reconstruction resolution. All these main effects together lead to a reduction factor of at most 0.2 in the statistical significance for pp → bbbbjj and of at most 0.5 for pp → bbγ γjj , as it can be seen in Fig. 18. The corresponding changes in the sensitivities to κ translate into the fact that, at this luminosity, the LHC will be sensitive to κ > 6.2 (7.7) at the 3σ (5σ ) level for the bbbbjj case and of κ > 7.7 (9.4) at the 3σ (5σ ) level for the bbγ γjj one. In the also interesting case of L = 3000 fb −1 , we get similar reduction factors. The reachable values of κ at this last HL-LHC stage via VBS configurations are predicted to be κ > 5.0 (6.3) at the 3σ (5σ ) level for the bbbbjj case and of κ > 6.1 (8.0) at the 3σ (5σ ) level for the bbγ γjj one.
In both cases we have seen that the values of κ that are closer to the SM value, say, κ ∈ [0.5, 2], are the most challenging ones to reach at the LHC. Even for the largest luminosity considered in this work, their corresponding statistical significances are always below 2σ . Hopefully, in this case gluon gluon fusion will be undoubtedly the only way to reach enough sensitivity, see refs. [17,21,22,25,36,38]. It is predicted, that, at L = 3000 fb −1 , statistical significances above 2σ will be always achieved in the gluon gluon fusion channel for λ ∼ λ SM .
The present study shows that double Higgs production via vector boson scattering is a viable and promising window to measure BSM deviations to the Higgs self-coupling and to deeply understand the scalar sector of the SM. Although all simulations are performed at the parton level, without hadronization or detector response simulation, and should be understood as a naive first approximation, we obtain very competitive results for the sensitivity to λ at the LHC. Because of this, we believe that the vector boson scattering H H production channel will lead to very interesting (and complementary to those of gluon-gluon fusion) findings about the true nature of the Higgs boson.