Probing the scalar potential via double Higgs boson production at hadron colliders

We present a sensitivity study on the cubic and quartic self couplings in double Higgs production via gluon fusion at hadron colliders. Considering the relevant operators in the Standard Model Effective Field Theory up to dimension eight, we calculate the dominant contributions up to two-loop level, where the first dependence on the quartic interaction appears. Our approach allows to study the independent variations of the two self couplings and to clearly identify the terms necessary to satisfy gauge invariance and to obtain UV-finite results order by order in perturbation theory. We focus on the $b \bar b \gamma \gamma$ signature for simplicity and provide the expected bounds for the cubic and quartic self couplings at the 14 TeV LHC with 3000 fb$^{-1}$ (HL-LHC) and for a future 100 TeV collider (FCC-100) with 30 ab$^{-1}$. We find that while the HL-LHC will provide very limited sensitivity on the quartic self coupling, precision measurements of double Higgs production at a FCC-100 will offer the opportunity to set competitive bounds. We show that combining information from double and triple Higgs production leads to significantly improved prospects for the determination of the quartic self coupling.


Introduction
Since the discovery in 2012 by the ATLAS and CMS collaboration [1,2], the Large Hadron Collider (LHC) has already disclosed an impressive amount of information on the properties of the resonance at 125 GeV, confirming so far the expectations of the Standard Model (SM).
The new particle is a narrow scalar [3,4], interacting with (third generation) fermions and vector bosons with a strength proportional to the mass of the particle [5,6]. All the expected main production and decay modes have been observed [1,2,[7][8][9][10][11][12]. Future runs at the LHC and future colliders will provide new information (such as the coupling to second generation fermions) and higher accuracy on the known couplings. The current measurements already indicate that New Physics (NP) effects cannot substantially affect the couplings of the Higgs boson with vector bosons and third generation fermions, placing the scale of NP well above the electroweak symmetry breaking scale.
The situation, however, is very different for the scalar potential on which we have not gained any relevant information so far and which is therefore largely unexplored. The reason is simply that the scalar potential, whose shape is ultimately responsible for Electroweak Symmetry Breaking (EWSB), can be probed only by measuring the Higgs self couplings. At hadron and lepton colliders, a direct sensitivity on the cubic or quartic Higgs self couplings can be achieved only via the simultaneous production of two or three Higgs bosons, respectively. Due to the smallness of the corresponding cross sections (around 32 fb in the case of pp → HH at 13 TeV [13][14][15][16] and 0.05 fb for pp → HHH [13,17]) these processes have not yet been observed at the LHC. Therefore the study of the Higgs self couplings is currently not only far from the precision level but also very challenging for the future.
In the case of double Higgs production, only exclusion limits are currently available and the most stringent result has been obtained by the ATLAS collaboration. Combining three different analyses (4b, bbτ τ, and bbγγ signatures) based on 27.5-36.1 fb −1 of data accumulated at 13 TeV [18][19][20][21], cross sections larger than 6.7 times the SM one can be excluded. This limit translates into the bound −5.0 λ SM 3 < λ 3 < 12.1 λ SM 3 , where λ 3 is the cubic coupling and λ SM 3 is its SM prediction. With a collected luminosity of 300 fb −1 , or even with 3000 fb −1 in the case of the High-Luminosity (HL) option, it is not still clear if the observation of SM production can be achieved. Although many phenomenological studies have been performed , the best experimental predictions for HL-LHC only provide upper limits on the SM cross sections.
In the case of λ 4 , the prospects are very uncertain. At the LHC, inferring information from triple Higgs production will be extremely challenging [45,46]. Its cross section is very small [13,17] and depends on the quartic interaction very weakly. Even a future 100 TeV proton-proton collider will need a considerable amount of integrated luminosity in order to obtain rather loose bounds [47][48][49].
Given the current and expected future results, new complementary strategies for the determination of the Higgs self couplings would be desirable. Recently, the possibility of probing the cubic Higgs self coupling λ 3 via precision measurements of single Higgs production channels at future lepton colliders [50] and at the LHC and future hadron colliders [51,52] has been suggested, exploiting the fact that next-to-leading order (NLO) EW corrections to the single Higgs production and decay modes involve λ 3 . The turning point for the possibility of determining the cubic interaction from single Higgs production measurements at the LHC has been the understanding that the different production channels depend on λ 3 in a very different way and that the effects are differential, the sensitivity being enhanced at threshold [52]. Even though the expected effects are small, a competitive sensitivity can be obtained by combining globally information from single Higgs measurements, total cross sections as well as distributions [52]. Since then, considerable effort has been invested in studying the feasibility of this strategy: predictions for the differential distributions for all the Higgs production channels have become available [53,54], and studies with more general (and realistic) scenarios for the existence of anomalous Higgs interactions [54][55][56][57] have appeared, also in combination with the direct double Higgs information [55][56][57][58]. Following the same logic, the λ 3 bounds have been extracted also from EW precision observables [59,60]. It is now clear that the indirect determination of λ 3 via precision measurements is expected to provide comparable bounds to those that are currently obtained via the direct searches for double Higgs production. Very recently it has been proposed that double Higgs production could be exploited for probing the quartic Higgs self coupling λ 4 via precise measurements [58,61]. The first studies at lepton colliders show that coarse bounds on λ 4 could be obtained and would complement the information from triple Higgs production, improving the ultimate results via a combination. At variance with the case of λ 3 , λ 4 -dependent loop corrections are ultraviolet (UV) divergent and in order to be renormalised they have to be performed in a Effective-Field-Theory (EFT) framework. The renormalisation procedure and the relevant counterterms have been provided in Ref. [58]. This framework has to be used also when the interest is focused only on independent variations of λ 3 and λ 4 , so that UV-finite results can be obtained.
The similar calculation for the case of hadronic collisions is computationally more involved, since the process pp → HH involves the loop-induced gg → HH partonic process at Born level and therefore the sensitivity on λ 4 originates from two-loop amplitudes. The first incomplete estimation of these effects has been presented in Ref. [62], showing the possibility of following this strategy also at future hadron colliders.
In this paper we analyse this strategy in detail and provide the first complete and consistent computation of the relevant contributions to gg → HH at two loops. All the two-loop diagrams involving λ 4 are taken into account and numerically evaluated without any further approximation via pySecDec [63,64]. Moreover, following the approach of Ref. [58], we take into account also corrections induced by additional λ 3 effects at two loops, which are non negligible for large values of λ 3 , and we renormalise the ensuing UV divergences. We perform this calculation at the differential level and we consider the bbγγ signature emerging from the decays of the Higgs bosons as a first application. This channel has been identified as the most promising one [26,42,[65][66][67][68][69] and it allows for the reconstruction of the di-Higgs invariance mass m(HH). Following the analyses in Ref. [26], we study the constraints that can be set on λ 3 and λ 4 via the measurement of the m(HH) distribution from bbγγ events for two different experimental setups: the LHC with 3000 fb −1 integrated luminosity (HL-LHC) and at a 100 TeV collider with 30 ab −1 integrated luminosity. The EFT parametrisation allows us to consider both the generic case, where λ 3 and λ 4 can vary independently, and a "well-behaved" EFT approach, where higher dimension operators induce smaller effects and λ 4 depends on λ 3 . In both cases we assume that the dominant BSM effects originate from the distortion of the Higgs potential, namely, anomalous interactions of the Higgs boson with other SM particles lead to subdominant effects. This approach is adequate to establish the sensitivity.
The paper is organised as follows. We first provide details on the computational framework, clarifying the theoretical assumptions, identifying the most relevant terms and describing the most important elements and features of the two-loop computation in sec. 2. Section 3 presents the results of the computation at the total as well as differential level, while in sec. 4 constraints that can be derived from future measurements at the LHC and at 100 TeV FCC are discussed in two different scenarios. We summarise our findings in sec. 5. Three appendices contain complementary and technical information.

Parametrisation of λ 3 and λ 4 effects
As already mentioned in the introduction, in order to vary the cubic and quartic Higgs self couplings λ 3 and λ 4 independently at all orders in perturbation theory in a consistent way, an EFT approach where operators are defined above the EWSB scale and respect all symmetries, hidden or not, has to be employed. This allows one to systematically identify gauge invariant and UV finite subsets of diagrams. For this reason, we will use the computational framework introduced and described in detail in Ref. [58]. In this section we summarise the most important aspects and we highlight some differences w.r.t. Ref. [58].
Starting from the SM Higgs potential we denote NP effects as V NP so that the general form of the potential can be written as where the symbol Φ refers to the Higgs doublet. Using the conventions of Ref. [70], the most general form of an SU (2)-invariant V NP potential reads One of the advantages of this parameterisation is that at tree-level λ 3 only depends on c 6 and λ 4 only on c 6 and c 8 . Indeed, after EWSB, we can rewrite V (Φ) as and thus define the self couplings λ 3 and λ 4 via are the values of λ 3 and λ 4 in the SM, respectively, and read In other words, the barred quantitiesc 6 andc 8 are simply c 6 and c 8 normalised in such a way that relations to κ 3 and κ 4 are simple. In particular  Using the parameterisation in eq. (2.3) and eqs. (2.5) and (2.6), or equivalently eqs. (2.8) and (2.9), we can trade κ 3 and κ 4 with only two other parameters,c 6 andc 8 . In so doing, we can always think of using the EFT approach as a way to obtain gauge invariant and UVfinite results in the anomalous coupling approach. 1 We note that, a priori, in a well-behaved EFT higher dimensional effects are expected to suppressed by a large scale Λ. Thus, in the first approximation, deviations in κ 3 and κ 4 are strongly correlated, i.e., (κ 4 −1) 6(κ 3 −1), see also eq. (2.9). Similarly to what as been done in Refs. [58,61,62], in this work we adopt as starting point an agnostic attitude towards the values that κ 3 and κ 4 can assume, in order to cover the sensitivity that future colliders can probe. We will later comment on bounds on κ 3 and κ 4 making different UV assumptions.
In this work we calculate the effects of anomalous cubic and quartic couplings in double Higgs production at hadron colliders. While λ 3 affects the gg → HH amplitude already at the Born level, λ 4 enters only via NLO EW corrections, i.e., at the two-loop level. Before discussing the details of the calculation it is convenient to anticipate what are the quantities that enter in our phenomenological predictions. In fig. 1 we display the one-loop diagrams of the Born amplitude in HH production. While the triangle (left diagram) depends on λ 3 , the box (right diagram) does not. Moreover, it is well known that the interference effects between the two diagrams leads to large cancellations. QCD corrections have been computed up to next-to-next-to-LO [14,16] and, besides reducing the scale dependence, they increase the LO cross section by roughly a factor of 2. In this work we will assume that QCD corrections factorise from the two-loop EW effects that we calculate. While the accuracy of this assumption has been directly tested only in very few cases [71][72][73], it has been often employed in the past, both due to the difficulty of calculating QCD-EW mixed corrections and due to the theoretical arguments supporting its validity. Two-loop corrections to HH production involve further λ 3 effects and introduce a λ 4 dependence, as can be seen in fig. 2. All the contributions arising from the two-loop topologies depicted in fig. 2 have been evaluated and renormalised via UV counterterms; more details concerning the calculation are given in Sec. 2.2.
Following the approach presented in Ref. [58] for e + e − collisions, we define the quantity 1 Note that using the alternative parameterisation V NP (Φ) ≡ ∞ n=3 c 2n Λ 2n−4 (Φ † Φ) n both λ3 and λ4 would depend on all the c i coefficients already at the tree level.  to be used in phenomenological investigations as where is the LO prediction. In eq. (2.11), σ 0 is the SM prediction, σ 1 corresponds to the leading contribution in the EFT expansion, being of order (v/Λ) 2 , while σ 2 is of order (v/Λ) 4 and corresponds to the squared EFT term. Clearly, no contribution proportional toc 8 appears at LO. The NLO corrections are included through the terms ∆σc 6 =c 2 6 σ 30c6 + σ 40c 2 6 +σ 20c 2 6 , (2.12) which are the loop corrections induced byc 6 on top ofc 6 and the two-loopc 8 -dependent part, respectively. They both originate from the topologies shown in Fig. 2. In the following we explain the rationale behind these formulae and the meaning of the different σ i(j)c i 6c j 8 terms entering them. First of all it is important to note that we organise the different contributions in terms ofc 6 andc 8 and not λ 3 and λ 4 . As explained in Ref. [58] this organisation reflects the necessary EFT expansion that has to be performed in order to renormalise UV divergences and obtain gauge invariant predictions. We recall thatc 6 can be directly related to λ 3 , whilē c 8 captures the violation of the relation κ 4 = 6κ 3 − 5, which holds if onlyc 6 is present, cf. eqs. (2.8) and (2.9).
Our goal is not to determine the ultimate precision that can be achieved at future colliders onc 6 andc 8 . Rather, we want to perform the first sensibility study on the determination of the cubic and quartic Higgs self couplings via double Higgs production at future hadron colliders. For this reason, SM EW corrections on top of σ LO are not taken into account. Since we are agnostic about the possible size ofc 6 , large cubic couplings are possible and lead to sizable enhancements via topologies such as (d) in fig. 2 [58,74]. For this reason, in ∆σc 6 we take into account all the contributions of orderc 3 6 andc 4 6 . These two contributions are relevant only for largec 6 , since otherwise they are suppressed w.r.t. the contributions appearing at LO. We remind that in Refs. [58,74] it has been shown that ∆σc 6 , and therefore σ pheno NLO , in general makes sense only in the range |c 6 | < 5. Outside this range perturbativity is violated for any prediction involving the bulk of HH production. We will comment more on this point in Sec. 3. At variance with Ref. [58], we include also the termσ 20c 2 6 in eq. (2.12). This term includes only part of the two-loop contributions of orderc 2 6 and its purpose is to preserve the large cancellations that are present in σ LO between different σ ic i 6 terms, since also in ∆σc 6 these cancellations are distributed among different σ i0 terms. On the other hand, it is relevant only forc 6 ∼ 2 where the cross section reaches the smallest value and the cancellations are the largest. 2 The quantity ∆σc 8 is the most relevant part of our computation and it solely induces the sensitivity onc 8 . At variance with Ref. [62], where only the topology (b) has been considered, in this term we take into account also all the contributions originating from topologies (e)-(i), which contribute at the same level and therefore cannot be ignored in any gauge-invariant calculation. 3 Also for the case ofc 8 , a theoretical bound based on the perturbativity requirement can be set [58] and corresponds to |c 8 | < 31.

Organisation of the calculation
In this section we give more details about our computational framework. Let us first consider the origin of the contributions in eqs. (2.12) and (2.13), in particular the presence of the termσ 20 . Using the same notations as for the σ i(j) terms, we define the different contributions of orderc i The SM term M 1L 0 receives contributions from both the one-loop triangle and box diagrams in fig. 1. The relation between eqs. (2.12) and (2.13) and the M i(j) terms is: In other words, ∆σc 6 and ∆σc 8 originate from the interference of M 1L with the terms with the largest dependence onc 6 , (c 2 6 M 2L 20 +c 3 6 M 2L 30 ), and all the terms that depend onc 8 , (c 8 M 2L 01 +c 6c8 M 2L 11 ). However, while the perturbative orders in ∆σc 8 and the interference terms emerging from the r.h.s. of eq. (2.17) are in one-to-one correspondence, this is not true for ∆σc 6 . The term 2 (M 1L 0 )(M 2L 20 ) * from the r.h.s. of eq. (2.16), which gives rise toσ 20 , multiplies the samec 6 powers as the term 2 (M 1L 1 )(M 2L 10 ) * , which we do not include in our computation. As already mentioned, we include the (formally subleading) termσ 20 because of the large cancellations among the triangle and box topologies at LO, and the fact that they contribute to differentc 6 powers; these cancellations are expected to be not substantially spoiled by NLO corrections. By keeping at the same level the entire M 1L amplitude of eq. (2.14) in the interference leading to ∆σc 6 we avoid that similar cancellations in NLO corrections are truncated by thec 6 expansion. As already mentioned, this is relevant only forc 6 ∼ 2, where σ has the minimum value, precisely due to the aforementioned cancellations. We remark, however, that this does not change the formal accuracy of our NLO corrections, which is of orderc 3 6 andc 4 6 . The two-loop contributions entering the different M ij sub-amplitudes can be further classified into three types: 3 Note that the topology (g) involves a H 5 interaction which in principle depends also on thec10 Wilson coefficient form the dimension-10 operator Φ † Φ − 1 2 v 2 5 . As discussed in Ref. [58], the effect of this diagram can be redefined as a constant shift onc6 and therefore our calculation is sensitive on a linear combination ofc6 andc10, which we set equal to zero. Nevertheless, thec6 andc8 contributions emerging from this diagrams are taken into account. See Ref. [58] for more details.
This classification is based on Feynman diagrams and can be easily understood from the topologies in fig. 2. The first category F corresponds to the factorisable topologies (d)-(j), together with the vertex counterterms in topologies (k) and (l). Their contributions are separately UV divergent, but their sum is finite, also for each separatec i 6c j 8 order considered in this work. In particular, the topologies (i), (j) and (l) can be evaluated together via the UV-finite P [HH] form factor given in Ref. [58], while all the remaining topologies from category F via the UV-finite V [HHH] form factor given in the same reference. We remind the reader that topology (d) is UV finite.
The non-factorisable two-loop contributions correspond to the topologies (a)-(c) which are not available. From a technical point of view, their computation is the most difficult and important part of this work. Details are given in sec. 2.3. Moreover, we find that numerically their phenomenological impact is non-negligible w.r.t. the factorisable ones.
We remind the reader that the Higgs wave-function renormalisation constant involves a quadratic dependence on κ 3 and therefore both a quadratic and linear dependence on c 6 [51,52]. Moreover, its contribution is UV-finite. Similarly to Ref. [58], its contribution is not included in the P [HH] and V [HHH] form factors and has to be separately added. The third category W corresponds to these additional contributions, which can be easily calculated via the LO diagrams and the SM contribution of λ 3 to the Higgs wave-function, namely, Based on the classifications we have just introduced, the different M 2L ij terms can be further divided into It will be useful to subdivide the amplitude further according to the spin exchanged in the s-channel. In view of the description of the calculation of two-loop non-factorisable diagrams, it is important to note that only the topology (a) includes both a spin-0 and spin-2 component; all the other topologies in fig. 2 are solely spin-0. In the case of one-loop diagrams, the triangle is also solely spin-0, while the box includes both a spin-0 and spin-2 component. Thus, the spin-2 contribution of the box diagram interferes only with the spin-2 component of the topology (a), while the spin-0 part of the box diagram and the triangle diagram interfere with all the two-loop topologies.
Since the diagrams in the topology (a), which involves both spin-0 and spin-2 components, lead to contributions of orderc 2 6 , we can further define where the first lower index denotes the spin component. With this notation we can directly express the M W 0,20 , M W 2,20 and also M F 30 terms as  In the topology (a) (double-box) there are in total 6 diagrams of which 3 are independent due to charge conjugation property of the fermion loop: 2 planar, (a 1 ) and (a 2 ), and 1 non-planar, (a 2 ). In the topology (b) (box-triangle) there are in total 3 diagrams, 2 of them are planar and charged conjugate, leading to (b 1 ), the remaining diagram is instead nonplanar, (b 2 ). The case of (c), is analogous to (b), including an H propagator. Topology (a), as we already said, contributes to both M N 0,20 and M N 2,20 , while topology (c) is in one-to-one correspondence with M N 10 . Topology (b) contributes also to M N 0,20 , which is therefore the only non-factorisable term receiving contributions from two different topologies. We can schematically summarise all this information as where the ⇐= arrow should be understood as "contributes to" and we have further remarked that M N 10 is all spin-zero. It is important to note that the sums of diagrams in each topology (a), (b) and (c) are separately finite and gauge invariant.  The calculation of all the non-factorisable two-loop diagrams is performed via numerical methods. As a first step, two-loop diagrams are generated with QGRAF [75] and the amplitudes are written in FORM [76] in d = 4 − 2 dimensions. Then, the amplitudes are projected onto spin-0 and spin-2 form factors. 4 Assigning the following on-shell p i momenta to the external particles, where all the p i are considered as incoming, both M 1L and M 2L , and any of their gaugeinvariant sub-amplitudes, can be projected onto two spin-0 and spin-2 bases A 1 and A 2 [77,78], and expressed via corresponding form factors denoted as F 0 and F 2 . Specifically, In eq. (2.28) 1 and 2 are the (transverse) polarisation vectors for the two incoming onshell gluons, while µ 1 and µ 2 (c 1 and c 2 ) are their corresponding Lorentz(colour) indices. The tensor bases A µ 1 µ 2 0 and A µ 1 µ 2 2 can be arbitrarily chosen and we decided to use the 4 In this work, this projection has been used also for the evaluation of F and W contributions.
orthonormal ones, which satisfy the relation 5 The two tensor bases A µ 1 µ 2 where p 2 T = (s 13 s 23 − m 4 H )/s 12 denotes the square of the Higgs-boson transverse momentum w.r.t. the gluons in the center-of-mass rest frame, using the convention s ij = (p i + p j ) 2 .
After the above projection, we obtain the spin-dependent non-factorisable amplitudes In other words, all the the N contributions entering our calculation can be expressed via the F N 0,20 , F N 2,20 and F N 0,01 form factors, which in turn depend on the non-vanishing spin-0 and spin-2 projections of the (a)-(c) topologies, F 0,a , F 0,b F 0,c and F 2,a .

Numerical evaluation of the form factors
The form factors F 0,a , F 0,b F 0,c and F 2,a are computed with pySecDec [63,64], a toolbox for the numerical evaluation of multi-loop integrals. We remind the reader that pySecDec can readily compute loop integrals with massive internal lines and/or off-shell legs. Moreover, compared to its predecessor SecDec 3 [79], it facilitates the creation of integral libraries, allowing for a direct incorporation of the code into the calculation of the full amplitude.
Before using pySecDec, we simplify the numerators of the the loop integrals in the form factors, in order to obtain tensor integrals that optimise the speed of the computation. It is important to note that the form factors F 0,a and F 2,a involve 7-propagator diagrams while F 0,b and F 0,c 6-propagator ones. Using propagator identities in FORM, we obtain a total of 11 integral expressions for F 0,a , 24 for F 2,a and 9 for F 0,b and F 0,c . The corresponding topologies are depicted in Appendix C. For simplicity, the overall coupling factors, colour factors (δ a 1 a 2 /2) and factor of (-1) due to fermion loop are removed from the tensor integrals. In particular, the quantities directly calculated via pySecDec areF 0,a ,F 2,a ,F 0,b andF 0,c , where we define The inner product stands for the contraction with polarisation vectors and summation over all physical polarisations. 6 The expression for the second projector in Ref. [77] contains a typo that is corrected here. F 2,a is related to F 2,a like in eq. (2.33) andF 0,c is related to F 0,c like in eq. (2.34).
In order to improve on the speed and convergence of the numerical evaluation, further measures are taken. First, only the finite parts are evaluated. To do this correctly, the integrands generated by pySecDec are multiplied with their prefactors, containing O( ) terms, before the integration. Nevertheless, we have cross-checked for specific phase-space points that UV-divergencies cancel for each diagram, although individual integral expressions can be separately UV-divergent.
Second, all integrals with the same denominator structure are added together before numerical integration. We have checked that the summation of several denominator structures prior to numerical integration does not lead to a faster convergence.
Third, different integrators were chosen for different integrals. A deterministic integrator like Cuhre [80], which is part of the Cuba library [81] and linked to pySecDec, is generally very fast and accurate for integrals with up to 5 integral dimensions. Beyond 5 dimensions, the integrator Vegas [82] is chosen. Furthermore, both Vegas and Cuhre give a χ 2 estimate, stating the probability that the uncertainty associated to the result is accurate. Tests have repeatedly shown that the Cuhre results can be trusted only if χ 2 is well below 1. Therefore, a routine was included to reperform the numerical integration with the more adaptive but generally slower integrator Vegas when the χ 2 value is too high. With this procedure we minimise cancellations and make sure that our numerical result is stable.
We have already mentioned that the UV finiteness of the form factors has been explicitly verified. Further tests have also been performed in order to ensure the correctness of the calculation. We have cross-checked the large m t limits for the (b) and (c) topologies (boxtriangle) against analytical results. By setting s 12 = m 2 H we have found perfect agreement with the expression given in Ref. [52]. Also, for the (a) topology (double-box), we have numerically tested that by artificially setting to m X the mass in the Higgs propagator connecting the two final-state Higgs, denoting the amplitude as M a,X , we obtain M a,X → M b [−2λ 2 3 /(λ 4 m 2 X )] in the limit m X → ∞. In other words, by integrating out the heavy state X, the (a) topology reduces to the (b) topology where as expected the quartic coupling is an effective coupling λ 4 = −2λ 2 3 /m 2 X . The factor of 2 originates from the number of diagrams contributing to the double-box amplitude, which is twice the number of diagrams contributing to the box-triangle amplitude.

Grids for phase-space integration
Up to this point we have discussed the strategy used for the evaluation of non-factorisable terms for a given phase space point. However, in this work we are interested in phenomenological predictions at colliders. Thus, the partonic squared matrix-elements have to be integrated over the phase-space and convoluted with parton-distribution-functions (PDFs). To this purpose, given the limited speed in the evaluation of the non-factorisable factors, it is helpful to build a grid that can be interpolated and quickly integrated over the relevant phase-space. In the following we explain how we have generated these grids, which have then been used together with an in-house Montecarlo for obtaining the phenomenological results of Sec. 4.  Let us start by discussing the spin-0 component at two loops. The box-triangle diagrams, topologies (b) and (c), depend on only one kinematic variable s 12 , hence a onedimensional grid for F 0,b and F 0,c is sufficient and with enough sampled values of s 12 a linear interpolation can be used. On the contrary, the double box diagrams, topology (a), depend on both s 12 and the angle θ between p 1 and p 3 . However, the dependence of F 0,a on θ is actually small and it can be approximated by the first few terms in the partial wave expansion [83] of F 0,a as We truncate the expansion in order to approximate the full results. We find that the θ dependence is weak, especially for s 12 < 4m 2 H , i.e., below the top-pair threshold in the loops. In this phase-space region the top-quark loop can be integrated out, obtaining an effective HHgg coupling among the Higgs bosons and the gluons. With such an EFT description in the m t → ∞ limit there is no θ dependence. Thus, the dominant contribution originates from the term without θ dependence, namely, the a 0 (s) term. In order to have the θ dependence under control and to test the validity of the partial wave expansion, we do not only include the first term but also the second term, 7 F 0,a (s, θ) ≈ a 0 (s) + a 2 (s)P 2 (cos θ) . (2.36) For each value of s, different values of θ have been sampled in order to perform a linear regression of a 0 (s) and a 2 (s). Afterwards, a linear interpolation is separately performed on both the values of a 0 (s) and a 2 (s). The validity of the truncation of the partial-wave expansion at a 2 (s) has also been investigated. First of all, we found that both the real and imaginary parts of a 2 (s) are substantially smaller than those of a 0 (s), as can be seen in fig. 4. Thus, contributions from higher-order a i (s) terms are expected to be even smaller than a 2 (s). Moreover we have estimated their contribution by comparing the value obtained with the approximation in eq. (2.36) after the regression and the actual value obtained. We can conclude that the truncation uncertainty is at the O(1%) level.  Table 2. Two-loop contributions to σ pheno NLO . We show for every entry the ratio with σ 0 at the same energy.
Let us conclude this section by commenting on the spin-2 contribution F 2,a . Although there is a large dependence on θ, we have verified that its contribution is strongly suppressed w.r.t. the spin-0 contribution. For this reason we safely ignore this contribution in our phenomenological study of Sec. 4.

Numerical Results
In this section we discuss the numerical results obtained for the m(HH) distribution and the total rates at different collider energies. The phenomenological analyses of Sec. 4 are based on these results.
In our calculation, we have used the following input parameters for the masses of the heavy SM particles, √ŝ , and we have used the Parton-Distribution-Functions (PDF) set CT14LO [84]. We remind the reader that in our calculation we renormalisec 6 in the MS scheme and we set the renormalisation scale to µ EFT = 2m H . Moreover, we assume both the Wilson coefficients c 6 andc 8 at the scale µ EFT . 0. 1. 1.

1.
2. In table 1, we list the three different σ i contributions entering the LO part of σ pheno NLO at 14, 27 and 100 TeV proton-proton collisions. Similarly, in table 2 we list all the twoloop σ ij contributions entering σ pheno NLO . We display in parentheses also their ratio with the LO prediction in the SM, σ 0 = σ SM LO . As can be seen in tables 1 and 2, cross sections considerably grow with the energy, while all the contributions induced byc 6 andc 8 mildly decrease in comparison with σ SM LO . Indeed, at large energies, the one-loop box diagrams is dominant w.r.t. the one with a triangle, which is the only one leading toc 6 contributions at LO and toc 8 contributions via loop corrections.
In fig. 5 we show four different contour plots for the 14 TeV energy. The upper plots show the ratio σ pheno NLO /σ SM LO , i.e., the ratio between our phenomenological prediction and the SM one, while the lower plots show the ratio σ pheno NLO /σ LO , which corresponds to the K-factor from two-loop corrections in our calculations. The left plots display these ratios in the (c 6 ,c 8 ) plane, while the right plots in the (κ 3 , κ 4 ) one. In the plots we consider the perturbativity regime |c 6 | < 5 and |c 8  total cross section. Forc 6 < 0 there is only a small dependence onc 8 , while forc 6 > 0 the dependence is sizable, and it even leads to negative cross sections for both large and positivec 6 andc 8 . These effects are induced by the loop corrections; the LO predictions cannot be negative since they originate from a squared amplitude. It can be seen also in the lower plots where the contour line for σ pheno NLO /σ LO = 0 is the same of σ pheno NLO /σ SM LO = 0 in the upper plots. For negative values our prediction is unphysical, so it cannot be used for phenomenological studies. This is caused by the sum ofc 6 andc 8 two-loop effects, which is large in absolute value. For the same reason also a region with σ pheno NLO /σ LO > 2 is present for large and positive(negative)c 6 (c 8 ). However, we do not exclude it since it is simply denoting a large one-loop K-factor.
We move now to the differential distributions. In fig. 6 we show the individual σ i (upper plot) and σ ij contributions (lower plots) to the m(HH) distribution at 14 TeV. 8 In the case of negative values we plot their absolute values and display the result as a dashed line. Moreover, we show in fig. 7 the ratio of any σ i and σ ij contribution over σ SM LO . In any plot this ratio is displayed as a black line, while we show in green the same result at the inclusive level, i.e., the values in parentheses in tables 1 and 2. We observe that thec 6 -andc 8 -induced contributions are most important close to threshold. Moreover, the quantities σ 1 , σ 30 and

General set up
In this section we discuss thec 6 andc 8 (κ 3 and κ 4 ) constraints that can be derived from the measurements of double Higgs production in proton-proton collisions at the LHC and a 100 TeV future collider. We consider the bbγγ signature, which has been identified as the most promising channel and allow for the reconstruction of the di-Higgs invariance mass m(HH).
In order to be close to a realistic experimental analysis, we follow the study of Ref. [26] for the case of HL-LHC and 100 TeV collisions with 30 ab −1 of luminosity. 9 We use the same selection cuts for the bbγγ signature, we divide the reconstructed m(HH) distribution in the same six bins and for each bin we take directly from Ref. [26] the predictions for the background and for the signal in the SM. Results in Ref. [26] take into account higher-order QCD corrections for both the signal and the background and also showering, hadronisation and detector effects. In our analyses we assume thatc 6 andc 8 effects factorise QCD corrections and we compute the effects of selection cuts (see Appendix A) adding H decays at the parton level. Thus, we also assume that showering, hadronisation and detector effects factorise the effect of selections cuts on the bbγγ signature.
In order to set limits onc 6 andc 8 we perform a χ 2 fit on the m(HH) distribution. For simplicity, as done in Ref. [26], we will include statistical uncertainties only. The impact of theoretical uncertainties and experimental systematic uncertainties is expected to be much smaller than statistical ones [26,62], therefore they would not in general lead to significant differences; some caveats are present for the 100 TeV case and will be discussed afterwards. On the other hand, we have found that assumingc 6 andc 8 effects as flat within each of the six bins of the reconstructed m(HH) distributions can strongly distort the results. Indeed, in each m(HH) bin,c 6 andc 8 effects are not flat over the full bbγγ phase-space. Thus, selection cuts have an impact not only on the total number of events observed but also on the ratio σ pheno NLO /σ SM LO . More details about the fit procedure can be found in Appendix B. Similarly to what has been in done in Ref. [58], we consider two different scenarios for setting bounds on Higgs self couplings: 1. Scenario 1: Well-behaved EFT (κ 3 = 1, κ 4 ∼ 6κ 3 − 1).
In Scenario 1, assuming that Nature corresponds toc 6 =c true 6 , we will analyse the constraints that can be set onc 6 . In the Scenario 2, settingc true 8 = 0, we explore the constraints that can be set on the (c 6 ,c 8 ) plane, for different value ofc true 6 . One may be tempted to study also a "Scenario 3", as done in Ref. [62], wherec 6 = 0 andc 8 = 0. However, this configuration is unstable. Indeed, it is easily spoiled by the running ofc 6 and c 8 at different scales, 10 since it is not protected by any symmetry and not emerging from an EFT expansion. For this reason we refrain from considering this scenario.

Scenario 1
We start considering the χ 2 function and the 1σ and 2σ bounds that can be obtained for c 6 assumingc true 6 = 0, at the HL-LHC and at a future 100 TeV collider. In fig. 9 we plot the χ 2 function 11 using σ pheno NLO or σ LO in the fit. Moreover, we show the relevance of fully differential information in the treatment ofc 6 andc 8 effects. In the case denoted as "flat µ-bin" in the plot, we assume that for each m(HH)-bin the impact ofc 6 effects can be evaluated via the ratio σ/σ LO without taking into account the selection cuts on the bbγγ final state, where σ can be either σ pheno NLO or σ LO . We remark that both in the "flat µ-bin" and normal cases, selection cuts are taken into account for the SM signal; the "flat µ-bin" concerns only the modelling ofc 6 andc 8 effects for the m(HH)-binning of the fit. More details are given in Appendix B. As can be seen in fig. 9, NLO effects, which in Scenario 1 corresponds to ∆σc 6 only, are relevant only for large values ofc 6 . On the contrary, the "flat µ-bin" assumption strongly distorts the χ 2 profile, especially for positive values ofc 6 . Indeed, as can be seen from the dashed lines, with this assumption the 2σ bounds at 14 TeV would be artificially improved. This effect is due to the fact that forc 6 2 the bulk of events is in the first bin(s) of the m(HH) distribution (see fig. 8) and selection cuts strongly depend on m(HH) especially close to the threshold (see Appendix A).

(4.2)
We now move to the case wherec true 6 can be different from zero. In fig. 10 we show 2σ bounds forc 6 as a function ofc true 6 . It turns out that ifc 6 is negative, bounds can be sizeably stronger. For instance, assumingc true 6 = −2 a limit −1.5. < κ 3 = 1 +c 6 < −0.5 can be obtained at HL-LHC, which is remarkably more stringent than in thec true 6 = 0 case of (4.1). In the case of 100 TeV, large and negative values ofc true under the "flat µ-bin" assumption as dashed lines. As can be seen, this assumption would have a strong effect to thec 6 bounds, especially forc 6 0.
In this context we want also to stress an important point that has been somehow overlooked in both theory and experimental studies on κ 3 -determination. In Fig. 11 we plot the 2σ constraints that can be obtained onc 6 by varying of σ exp /σ SM , where σ exp is the measured value and σ SM is the SM prediction. We derive the constraints using two different approximations: σ pheno NLO and σ LO . As can be seen, for |c 6 | 5, where perturbativity is violated, the constraints onc 6 strongly depend on the choice between σ pheno NLO and σ LO . When data are fitted with σ LO predictions,c 6 or equivalently κ 3 is a parameter of ignorance that only for |κ 3 − 1| = |c 6 | 5 coincides to the quantity one is interested in. Outside this range,c 6 (or κ 3 ) is only suggesting how far from the SM predictions is the experimental result. The usage of σ pheno NLO or any higher-order corrections in the place of σ LO is not improving this situation, since the regime is not perturbative for |c 6 | 5. In conclusion, one can set bounds outside the |κ 3 − 1| = |c 6 | 5 range, but only within this region they properly refer to the quantities we are interested in and defined via parameters in the Lagrangian.

Scenario 2
This scenario allows us to discuss the most important phenomenological results of this work, i.e., the expected constraints onc 6 andc 8 (κ 3 and κ 4 ) that can be obtained via double Higgs production at HL-LHC and a 100 TeV future collider. Assumingc true 8 = 0, these constraints are shown in the left and right plot of Fig. 12, respectively. We show 2σ results and again the effect due to the "flat µ-bin" assumption, the red area corresponds to the region where the cross-section is negative (cf. left plots in Fig. 5). As already mentioned, no phenomenological study can be performed in this configuration. Similarly, for a given (c 6 ,c 8 ), predictions for some bins can be negative, while positive for others; we retain the information only for those bins where the cross-section is predicted to be positive. As can be seen from fig. 12, at HL-LHC the presence ofc 8 contributions is not sizeably affecting the result in (4.1), obtained under the assumptionc 8 = 0. On the other hand, no sensible constraints can be obtained at the HL-LHC on thec 8 parameter. In other words, with a complete calculation and taking into account selection cuts and background effects, we find a much less optimistic result than in Ref. [62].
Results at 100 TeV collisions are qualitatively very different than at the HL-LHC. The bounds onc 6 are affected by the presence onc 8 . As can be seen from the right plot of Fig. 12, the bounds are 0.4 < κ 3 = 1 +c 6 < 2, which is less precise than (4.2), obtained under the assumptionc 8  lead to strong constraints in the (c 6 ,c 8 ) plane. However, we remind the reader that we do not take into account theory and experimental systematic uncertainties. As said for the corresponding results in Scenario 1, these results may be affected by the aforementioned uncertainties.
Last but not least, in Fig. 15 we compare the constraints obtained forc true 6 = 0 at 100 TeV (right plot of Fig. 12) with the corresponding ones obtained following the analysis in Ref. [65,85], based on the bbbbγγ signature emerging from pp → HHH production. 12 Triple Higgs bounds are derived via two different assumptions on b-tagging efficiency: optimistic (80%) and conservative (60%). As can be seen in Fig. 15, double Higgs bounds are stronger than those from triple Higgs with the optimistic assumption. Especially, they are complementary to those from triple Higgs with the conservative assumption and their combination can lead to stronger results. We also show the corresponding comparison in (κ 3 , κ 4 ) plane taking into account the perturbative bounds onc 6 andc 8 .   Figure 15. Comparison between expected 2σ-bounds from HH(bbγγ) and HHH(bbbbγγ) at 100 TeV. The right plot in (κ 3 , κ 4 ) plane takes in account the perturbativity bounds onc 6 andc 8 .

Conclusion
The experimental determination of the Higgs potential and in particular of the Higgs self couplings is one the most far fetching goals of the HL-LHC and of future colliders. Its importance is matched only by the difficulty of such an endeavour: rates for multiple Higgs production which are directly sensitive to the self couplings, are very low making it hard to study distributions where most of the sensitivity actually lies. This is certainly true for the cubic coupling at the LHC, which can be accessed directly via HH production, but becomes dramatic for the quartic coupling: its direct determination calls for measurements in the HHH final state, whose production rate will be small even at a future 100 TeV pp colliders.
The challenge on the one hand and the high-stakes on the other hand have provided strong motivation to the theoretical and experimental high-energy-physics community to devise alternative strategies. Among them, a new approach has emerged building up from the simple idea that single Higgs cross sections might display a sensitivity on the cubic coupling at higher orders. Since the first proposal in the context of future e + e − colliders [50], the idea has been developed and extended to hadron colliders, eventually proving to be competitive with the direct determinations. A very first experimental analysis by CMS [86] based on the proposal of refs. [52,54] has confirmed the expectations of the theoretical studies.
Recently, some of us have proposed to extend the idea further and determine the (cubic and) quartic coupling exploiting the sensitivity coming from loop effects in HH in the context of future e + e − colliders [58] . In this work we have moved one step further and explored the reach of hadron colliders by determining the sensitivity to the (cubic and) quartic coupling of the main double Higgs production channel, gg → HH, up to two loops. Being a technically challenging two-loop computation we have employed the most up-todate numerical multi-loop techniques, providing for the first time a complete and consistent calculation of these effects.
We have considered two different scenarios, one "EFT-like" where the cubic and quartic couplings are related and one where they are varied independently. Our results clearly indicate that while the HL-LHC will have limited sensitivity, at the FCC-100 the precision on HH differential measurements will be such that HH will be more sensitive to independent deviations in the self couplings than HHH production itself. The best constraints on the quartic will therefore be obtained by combining HH precision measurements with the direct determinations from HHH.

MOVE-IN Louvain
Cofund grant. The work of X.Z. is supported the European Union's Horizon 2020 research and innovation programme as part of the Marie Skłodowska-Curie Innovative Training Network MCnetITN3 (grant agreement no. 722104). This work has received funding from the ERC grant "MathAm" and from F.R.S.-FNRS under the 'Excellence of Science' EOS be.h project n. 30820817. Computational resources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI).

A Cut efficiency
In this section we explicitly write the cuts used in our analysis. The cuts are the same as in Ref. [26], on which our analysis is based. Specifically, at 14 TeV, we have, In fig. 16, we show the differential cut efficiency for the signal, assuming SM double Higgs production and narrow-width approximation. In other words we plot the ratio between the number of events predicted in the SM with and without the cuts as function of m(HH). Since spin-0 contributions dominate for both SM and BSM cases, cut efficiencies for BSM cases are very similar.
The zero efficiency in the 250 GeV < m(HH) < 300 GeV phase-space region is not a surprise; when Higgs boson pairs are produced at the threshold, both the bb and γγ pairs from the Higgs decays are back-to-back and therefore rejected by the cuts ∆R(b, b) < 2 and ∆R(γ, γ) < 2. Increasing the energy, both Higgs can have non-vanishing transverse momentum and therefore their decay products can be not back-to-back and tend to be collimated for very high energies.

B Fit details
In this Appendix we describe in detail the χ 2 functions that have been used in this work for extracting from m(HH) distributions 1σ and 2σ bounds in the (c 6 ,c 8 ) parameter space.
The general formula of the χ 2 that has been exploited for our results has two degrees of freedom,c 6 andc 8 , and reads Bounds onc 6 andc 8 have been obtained following a fit procedure similar the one presented in Ref. [26], from which we have taken also the selection cuts (see Appendix A) and binning in the m(HH) distribution. For this reason, the value of N BKG i is directly taken from Ref. [26]. On the contrary, N HH i (c 6 ,c 8 ) is derived from the value N HH i (0, 0), the SM prediction, from the same reference, which takes into account also higher-order QCD corrections. Assuming that these effect factorise withc 6  where in the right equation we have understood the dependence onc 6 andc 8 and σ SM LO = σ pheno NLO |c 6 =0,c 8 =0 . The quantity Φ i corresponds to the bbγγ phase space such that the reconstructed m(HH) value is within the bin i. Within all the work, unless differently specified, we take into account the selection cuts of Appendix A in Φ i . When we say "flat µ-bin" we precisely refer to the case where selection cuts are not taken into account for the definition of µ theory i .

C Topologies of the integral expressions from non-factorisable two-loop contributions
In this Appendix we show the topologies of the integral expressions obtained from nonfactorisable two-loop contributions. Each topology can lead to more than one integral expression. In Fig. 17 we show those relevant forF 0,a andF 2,a , while in Fig. 18 those forF 0,b , which are relevant also forF 0,c . Thick lines correspond to massive propagators, dashed with mass m H while solid with mass m t . The solid-thin lines corresponds to massless propagators.