Effective field theory analysis of double Higgs boson production via gluon fusion

We perform a detailed study of double Higgs production via gluon fusion in the effective field theory (EFT) framework where effects from new physics (NP) are parametrized by local operators. Our analysis provides a perspective broader than the one followed in most of the previous analyses, where this process was merely considered as a way to extract the Higgs trilinear coupling. We focus on the hh → bb̄γγ channel and perform a thorough simulation of signal and background at the 14 TeV LHC and a future 100 TeV proton-proton collider. We make use of invariant mass distributions to enhance the sensitivity on the EFT coefficients and give a first assessment of the impact of jet substructure techniques on the results. The range of validity of the EFT description is estimated, as required to consistently exploit the high-energy range of distributions, pointing out the potential relevance of dimension-8 operators. Our analysis contains a few important improvements over previous studies and identifies some inaccuracies there appearing in connection with the estimate of signal and background rates. The estimated precision on the Higgs trilinear coupling that follows from our results is less optimistic than previously claimed in the literature. We find that a ∼30% accuracy can be reached on the trilinear coupling at a future 100 TeV collider with 3 ab−1. Only an Oð1Þ determination instead seems possible at the LHC with the same amount of integrated luminosity.


I. INTRODUCTION
Unveiling the dynamics at the origin of electroweak symmetry breaking (EWSB) is a primary goal of the experiments performed at the Large Hadron Collider (LHC) at CERN.The discovery at Run1 of a new boson with mass m h ≃ 125 GeV and properties similar to those predicted for the Standard Model (SM) Higgs boson has been a major leap forward in this direction [1,2].The many direct searches and precision measurements performed by the ATLAS and CMS Collaborations all indicate the existence of a gap between the electroweak (EW) scale and the scale of NP, unless the latter is very weakly coupled to the SM sector.This justifies the use of effective field theory (EFT) to give a low-energy parametrization of NP effects in terms of a series of local operators.Measuring the coefficients of these effective operators would give access to a wealth of information on the UV dynamics, most importantly whether it is strongly or weakly coupled.Higgs studies at Run1 have mostly focused on on-shell single-production and decay processes, thus probing the strength of the underlying EWSB dynamics at the scale Q ¼ m h .They set limits on possible modifications of the Higgs couplings, whose naive expected size is of order δc=c ∼ ðg 2 Ã =g 2 SM Þm 2 h =m 2 Ã , where m Ã is the mass of the new states, g Ã is their coupling strength to the Higgs boson, and g SM is a SM coupling.The improved performances in energy and luminosity that will characterize the LHC Run2, on the other hand, give the opportunity to directly probe the EWSB dynamics at much higher energies (Q ∼ E ≫ m h ) through the study of 2 → 2 scattering processes.For a typical scattering energy E, effects from NP are expected to be of order ∼ðg 2 Ã =g 2 SM ÞE 2 =m 2 Ã and hence enhanced by a factor E 2 =m 2 h compared to those entering on-shell Higgs processes.Exploring higher energies thus gives access to potentially larger corrections, but at the same time poses the issue of assessing the validity of the EFT description.Determining at which point this latter breaks down, in fact, requires adopting a power counting to estimate the size of the local operators in terms of the parameters (masses and couplings) characterizing the UV dynamics.At the same time, the power counting puts the limits on the effective operators into perspective and helps infer how much theoretical space is being probed.Double Higgs production via gluon fusion is one example of a scattering process that can disclose key information on the EWSB dynamics, its underlying symmetries, and its strength.It is the only process potentially observable at the LHC that can give access to the quartic couplings among two Higgs bosons and a pair of gluons or of top quarks, as well as to the Higgs trilinear self-coupling.Previous studies of this reaction mostly focused on the extraction of the trilinear coupling in the context of the SM [3][4][5][6][7][8][9][10][11].An analysis based on the wider EFT perspective, however, can give access to a much richer spectrum of information on the UV dynamics.In this work we provide such analysis at the 14 TeV LHC and a future 100 TeV proton-proton collider by focusing on the hh → b bγγ final state.We perform a detailed Monte Carlo (MC) simulation of signal and background, and we use the kinematic distributions to maximize the sensitivity on different effective operators.Previous studies in the EFT context appeared in Refs.[12,13], and a first analysis of the impact of NP on the kinematic distributions can be found in Ref. [14].
The paper is organized as follows.Section II contains a detailed discussion of the parametrization of double Higgs production within the framework of EFT.The power counting of the relevant coefficients is reviewed and a few scenarios are examined where NP gives large modifications to the Higgs trilinear coupling, whereas other couplings stay close to their SM values.The relevance of dimension-8 operators is also investigated, and the validity of the EFT description is assessed.Section III reviews the phenomenology of double Higgs production at the LHC and a future 100 TeV collider, giving a first estimate of the sensitivity on the kinematic tail at large invariant masses of the Higgs pair.Our Monte Carlo analysis of the b bγγ decay mode is illustrated in Sec.IV, while results on the sensitivity to the coefficients of the effective Lagrangian are collected in Sec.V. We compare with previous analyses of the b bγγ final state in Sec.VI.Section VII puts our study into context by briefly discussing existing studies of other final states and gives an outlook on possible future developments.We collect useful formulas and further details on our simulation of the b bγγ background into the Appendixes.

II. EFFECTIVE PARAMETRIZATION AND POWER COUNTING
Corrections due to the exchange of new heavy states can be conveniently parametrized by means of low-energy effective Lagrangians.There are two formulations that are suited to the study of Higgs physics.The first one assumes that the Higgs boson is part of a weak doublet, as in the SM, and that SUð2Þ L × Uð1Þ Y is linearly realized at high energies.It thus is referred to as the "linear" Lagrangian.In the second, more general formulation, SUð2Þ L × Uð1Þ Y is nonlinearly realized, hence the name of "nonlinear" Lagrangian, and the physical Higgs boson is a singlet of the custodial symmetry, not necessarily part of a weak doublet.Both parametrizations have been reviewed in Ref. [15], to which we refer the reader for a more in-depth description.The experimental data collected at the LHC during Run1 seem to indicate that the couplings of the newly discovered boson have values close to those predicted for the SM Higgs.Although such experimental information is still preliminary and awaits the confirmation of Run2, it motivates the use of the linear Lagrangian for future studies.Indeed, small deviations from the SM can be obtained if the Higgs boson belongs to a doublet, provided the new states are much heavier than the weak scale.The nonlinear formulation is still useful, however, in those cases where large deviations in the Higgs couplings are allowed.As we will see later in our analysis, this is especially true for double Higgs production, where additional couplings not accessible via single Higgs processes can be extracted.
In the linear Lagrangian, the operators can be organized according to their dimension, The lowest-order terms coincide with the usual SM Lagrangian L SM , whereas L n contains the deformations due to operators of dimension n, with n > 4. 1 For our purposes it is sufficient to focus on the operators involving the Higgs boson.The relevant ones in L 6 are 2 GeV and m h ¼ 125 GeV is the physical Higgs mass.To classify the various operators in our effective Lagrangian and estimate their coefficients, we adopt the Strongly Interacting Light Higgs (SILH) power counting [16].This is based on the assumption that the NP dynamics can be broadly characterized by one mass scale m Ã , at which new states appear, and by one coupling strength g Ã .The latter, in particular, describes the interactions between the Higgs boson and the new states.When building higher-order operators starting from the SM Lagrangian, each extra insertion of the Higgs doublet is weighted by a factor 1=f ≡ g Ã =m Ã , while each additional covariant derivative is suppressed by m Ã .If the Higgs is a composite Nambu-Goldstone (NG) boson [17], the operators O 6 , O u , and O g can be generated only through some small explicit breaking of the global invariance since they 1 Here and in the following we assume that baryon and lepton numbers are conserved by the NP dynamics and thus by the series of higher-dimensional operators.The only dimension-5 operator invariant under the SM gauge symmetry violates the lepton number and can be omitted.
violate the Higgs shift symmetry [16]. 3In this case the SILH estimate of the coefficients in Eq. ( 2) is where λ is a weak spurion coupling suppressing cg , while the suppression of c6 and cu is controlled, respectively, by ðm 2 h =2v 2 Þ and y t .In realistic models, the minimum amount of explicit breaking entering cg is given by the top Yukawa coupling y t , so that one expects y t ≲ λ ≲ g Ã .For example, one has λ ¼ y t in models with fully composite t R , whereas partial compositeness of both t L and t R leads to λ ¼ ffiffiffiffiffiffiffiffi g Ã y t p [16].In the case of the nonlinear Lagrangian, the terms relevant for our analysis are where h denotes the physical Higgs field (defined to have vanishing vacuum expectation value).Compared to the linear Lagrangian, the couplings c i of Eq. ( 4) effectively resum all the corrections of order ðv 2 =f 2 Þ.Therefore, the nonlinear Lagrangian only relies on the derivative expansion, in which higher-order terms are suppressed by additional powers of ðE 2 =m 2 Ã Þ.The linear formulation further requires ðv=fÞ < 1 for the expansion in powers of the Higgs doublet to be under control.When both parametrizations are valid, the coefficients of the two Lagrangians are related by the following simple formulas: where α 2 ≡ g 2 =4π.Notice in particular that the same operator O u that gives a (nonuniversal) modification of the top Yukawa coupling also generates the new quartic interaction tthh.
It is worthwhile at this point to comment about the normalization of the operator O 6 in Eq. ( 2).This differs from the one of Ref. [15], where it was defined ΔL ⊃ −c 0 6 λ 4 ðH † HÞ 3 =v 2 [here we introduce the symbol c0 6 to distinguish it from c6 in Eq. ( 2)], with λ 4 denoting the coefficient of the ðH † HÞ 2 operator in the SM Lagrangian.We find , which should be compared with the relation between c 3 and c6 (valid at all orders in c6 ) appearing in Eq. (5).For small values of c0 , so that the dependence of c 3 upon c0 6 is the same as in Eq. ( 5) up to small corrections.The parametrization of Eq. ( 2) is thus more convenient than that of Ref. [15] when c6 (or c0 6 ) becomes large, since the formula for the trilinear coupling c 3 remains linear in c6 .This is relevant for this analysis, since we anticipate that the results of Sec.V will constrain values of c6 larger than 1 at the LHC.

A. Modified power counting for the Higgs trilinear coupling
The estimates of Eq. ( 3) show that assuming the SILH power counting the modification of the Higgs trilinear coupling is expected to be of order ðv=fÞ 2 , i.e., of the same size as the shifts to other Higgs couplings.These latter, however, are already constrained from single-Higgs measurements to be close to their SM value, in particular, the coupling of the Higgs to two vector bosons.It is thus interesting to ask whether there are scenarios, characterized by a power counting different from the SILH one, where the Higgs trilinear coupling can get a large modification from NP effects while all the other Higgs couplings are close to their SM values, in agreement with the current LHC limits.Interestingly the answer to this question is affirmative: it is possible to imagine at least a few scenarios where the largest NP effects are in the trilinear coupling.
A first possibility is a scenario where the Higgs is a generic bound state-i.e., not a Nambu-Goldstone bosonof some new strong dynamics.In this case no weak spurion suppression is required to generate O 6 , and the naive estimate of c6 is enhanced by a factor g 2 Ã = λ4 compared to the SILH case, where we conveniently defined λ4 ≡ m 2 h =2v 2 ≃ 0.13.The Higgs trilinear coupling (in SM units) is thus expected to be of order c 3 ¼ 1þ Oðg 2 Ã v 4 =f 2 m 2 h Þ, and large modifications are possible for g Ã large even if ðv=fÞ 2 ≪ 1.For example, ðv=fÞ 2 ¼ 0.05 gives cH ∼ 0.05 and c6 ∼ 3.5ðg Ã =3Þ 2 , corresponding to c 3 ≃ 0.93 þ 3.5ðg Ã =3Þ 2 .The price to pay for this enhancement, however, is the tuning, which is now required to keep the Higgs mass light, since naively one would expect m 2 h ∼ m 2 Ã .This is in addition to the Oðv 2 =f 2 Þ tuning that must occur in the vacuum alignment even for a NG boson Higgs.Notice that similar to O 6 , other operators that previously required a breaking of the Goldstone symmetry to be generated will now be unsuppressed.In particular, cg and the coefficient of the operator H † HB μν B μν get enhanced by a factor ðg 2 Ã =λ 2 Þ compared to their SILH estimate, thus leading to modifications of the Higgs couplings to gluons and photons of order v 2 =f 2 times their (loop-induced) SM value.However, for ðv=fÞ 2 ≪ 1 the current LHC constraints on these couplings can easily be satisfied.
It is possible to avoid the tuning of the Higgs mass while keeping an enhanced trilinear coupling by envisaging a different scenario.Consider a theory where a new strongly interacting sector, characterized by a mass scale m Ã and a coupling strength g Ã , couples to the SM sector only through the Higgs mass term portal: L int ¼ λH † HO.The Higgs field is thus elementary and couples to the composite operator O, made of fields of the strong sector, with coupling λ.This leads to the following estimates: where Δλ 4 is the correction to the coefficient of the operator ðH † HÞ 2 induced by the strong dynamics.The contribution to cu follows from a Higgs-dependent modification of the top quark kinetic term arising at the 1-loop level, and it is subleading compared to cH .Similar effects are also generated for the lighter quarks proportional to their Yukawa coupling squared.Corrections to Eq. ( 6) of higher order in λ are suppressed by powers of ðλv 2 =m 2 Ã Þ.An additional factor g 2 Ã =16π 2 should be included in the above estimates if they arise from diagrams involving a loop of strongly coupled particles. 4Since c6 ∼ ðλ= λ4 Þc H , one can have an enhanced NP correction to the Higgs trilinear coupling for λ > λ4 .However, requiring that no fine-tuning occurs in the coefficient of the quartic operator-i.e., in the physical Higgs mass-implies λ ≲ ffiffiffiffi ffi λ4 p g Ã .Furthermore, only for ðλv 2 =m 2 Ã Þ ≲ 1 one can make sense of the series in powers of the Higgs field, since when this inequality is saturated all orders in λ become equally important and perturbation theory is lost.From the above considerations it follows the bound c6 ≲ Oð1Þ.In this scenario it is thus possible to avoid tuning the Higgs mass, but the corrections to the Higgs trilinear are at most of Oð1Þ; larger enhancements require fine-tuning.For example, g Ã ¼ 4π, m Ã ¼ 500 GeV, and λ ¼ 4π ffiffiffiffi ffi λ4 p ≃ 4.5 give cH ∼ 0.03 and c6 ∼ 1.These examples show that it is possible, in scenarios characterized by a power counting different from the SILH one, to have large modifications to the trilinear coupling while keeping the other Higgs couplings close to their SM values.In some cases, like that of a generic composite Higgs, c6 can be of order 1 or larger without invalidating the perturbative expansion in powers of the Higgs field.

B. The relevance of higher-order operators
Let us now investigate more quantitatively the validity of our effective Lagrangians and discuss the importance of higher-order operators.We will do so by adopting the SILH power counting.We already mentioned that in the case of the linear Lagrangian, for a given process, operators with higher powers of the Higgs doublet imply corrections suppressed by additional factors ðv=fÞ.This means that if v=f ∼ 1, the linear Lagrangian cannot be consistently used anymore, although we can still rely on the nonlinear one.As regards the derivative expansion, the same considerations apply instead to both descriptions.The perturbative parameter in that case is ðE=m Ã Þ, so that naively higher-order operators become important at E ∼ m Ã , where the effective description breaks down.However, in certain processes it can happen that the leading operators are suppressed due to symmetry reasons, and the higher-order ones become relevant at large energies below the cutoff scale.This is, in fact, the case of double Higgs production via gluon fusion, where the dimension-6 operator O g is suppressed by two powers of a weak spurion if the Higgs is a NG boson.
The dimension-8 operators relevant for double Higgs production via gluon fusion are They mediate scatterings of two gluons into two Higgs bosons with total angular momentum equal to, respectively, 0 and 2, and can therefore lead to different kinematic configurations, as we will discuss later on.Both operators, however, are similar from the power-counting viewpoint.
Neither of the two breaks the shift symmetry of a NGB Higgs, and thus no spurion suppression appears in the estimate of their coefficients, We can now compare the contributions of the various operators to double Higgs production via gluon fusion. 6By using the estimates of Eqs. ( 3) and ( 8), the scattering amplitude can be schematically expressed as 4 This is what happens, for example, when the strong dynamics consists of a single real scalar field S, neutral under the SM gauge group.By requiring invariance under the parity S → −S, the Lagrangian describing the strong dynamics is where The first term in the square parentheses of Eq. ( 9) corresponds to the SM contribution plus the Oðv 2 =f 2 Þ correction implied by O H and O u . 7The second and third terms in the parentheses denote the contributions of, respectively, O g and O gD0 , O gD2 , and thus define the coupling strengths of the interactions mediated by these operators.Since these couplings grow with E, at sufficiently large energies, yet below the cutoff scale, they can become larger than the top Yukawa coupling.It is thus possible to obtain NP corrections larger than the SM contribution within the validity of the effective Lagrangian.This should be contrasted with the corrections from dimension-6 operators to the on-shell production and decay rates of the Higgs boson: those are at most of order Oðv 2 =f 2 Þ, hence always smaller than the SM term.We thus see the potential advantage of studying 2 → 2 scatterings at high energies compared to the single-Higgs measurements performed during Run1 at the LHC: going off-shell at higher energies one can directly probe the strength of the NP dynamics in the regime in which it gives large effects [18].In practice, as we will discuss more in detail in the following, such a regime is difficult to exploit in the process gg → hh → b bγγ considered in this work.This is due to the steep falloff of the gluon Parton Distribution Function (PDF) at large x and to the small final-state branching ratio, which strongly limit the exploration of the region with large hh invariant mass.For λ > ffiffiffiffiffiffiffiffi g Ã y t p the coupling g 6 ðEÞ is the first to become larger than y t at a scale E ∼ y t fðg Ã =λÞ.The smaller λ is, the higher the crossover scale becomes, until E ∼ y t fðg Ã =y t Þ 1=2 is reached where g 8 ðEÞ ∼ y t .Above this energy, for λ ≲ ffiffiffiffiffiffiffiffi g Ã y t p , the largest effects come from the dimension-8 operators.The validity of the effective Lagrangian extends up to energies E ∼ m Ã , so that the coupling strengths are bounded by As we anticipated, studying the kinematic region where higher-order operators give a contribution larger than the SM one is challenging in double Higgs production via gluon fusion.However, it is still interesting, and relevant for the analysis carried out in this work, to ask at which scale the dimension-8 operators become more important than the dimension-6 ones, independently of their absolute size.From Eq. ( 10) it follows that g 6 ðEÞ ∼ g 8 ðEÞ for E ∼ λf; at this scale g 6 ðEÞ ∼ λ 2 =g Ã .Hence, a further condition to be satisfied in an analysis that includes only dimension-6 operators neglecting those with dimension 8 is Let us then indicate with δ the precision obtained in such an analysis on c 2g ¼ cg ð4π=α 2 Þ, by making use of events with invariant masses up to Ē.This means that the smallest value of g 6 probed is g min ∼ ffiffi ffi δ p ð Ē=vÞ, so that Eqs. ( 11) and ( 12) can be recast into the following constraints: These inequalities define the region in the ðg Ã ; λÞ parameter space, sketched in Fig. 1, which can be sensibly probed, within the validity of the effective theory, by an analysis including only dimension-6 operators.Clearly, smaller values of g min , obtained by increasing the precision δ for a fixed energy Ē, lead to a larger viable region.Notice that the two cases λ ¼ y t (fully composite t R ) and λ ¼ ffiffiffiffiffiffiffiffi g Ã y t p (partially composite t L and t R ) can be probed only if g min < y t .This is not the case, of course, if the analysis does not have enough sensitivity to probe the SM cross section.
FIG. 1 (color online).Region in the plane ðg Ã ; λ=g Ã Þ, defined by Eqs. ( 13) and (14), that can be probed by an analysis including only dimension-6 operators (in white).No sensible effective field theory description is possible in the gray area (λ < g min ), while exploration of the light blue region (g min < λ < ffiffiffiffiffiffiffiffiffiffiffiffi g Ã g min p ) requires including the dimension-8 operators.

C. Cross section of double Higgs production
We can now discuss our parametrization of the cross section of double Higgs production via gluon fusion.We will use the nonlinear Lagrangian (4) and start by neglecting higher-derivative terms (which correspond to dimension-8 operators in the limit of linearly realized EW symmetry).The effect of the neglected derivative operators will then be studied by analyzing their impact on angular differential distributions and shown to be small in our case due to the limited sensitivity on the high m hh region.
The Feynman diagrams that contribute to the gg → hh process are shown in Fig. 2. Each diagram is characterized by a different scaling at large energies ffiffi ffi ŝ p ¼ m hh ≫ m t ; m h .We find where A □ , A △ are the amplitudes of, respectively, the box and triangle diagram with Higgs exchange, A △nl is the amplitude of the triangle diagram with the t thh interaction, and A 3 and A 4 are the amplitudes of the diagrams with the Higgs-gluon contact interactions.One can see that each NP contribution affects the m hh distribution in a different way.In particular, the diagrams that depend on the Higgs trilinear coupling c 3 are always suppressed at large ŝ, and their contribution affects the process mostly at threshold.Modified values of the top Yukawa coupling c t and the nonlinear interactions c 2t and c 2g , instead, tend to increase the cross section at higher invariant masses.Finally, including the dimension-8 operators would lead to an additional contribution to A 4 growing as ŝ2 and distort the tail of the m hh distribution.A shape analysis can thus help to differentiate the different effects and break the degeneracy of the total cross section on the Higgs couplings.This will be our strategy in the study of double Higgs production discussed in the next section, where we will use m hh as the main kinematic variable to characterize signal events.By focusing on gluon fusion, the total cross section for the process pp → hh can be written as a simple polynomial of the parameters of the effective Lagrangian, 8 The LO value of the numerical coefficients A i and of the SM cross section σ SM is reported in The contribution from vector boson fusion is smaller by at most 1 order of magnitude and can be isolated by selecting the number of jets in the final state.See, for example, Refs.[19,20] for up-to-date studies with m h ¼ 125 GeV. 9 We computed the A i 's by performing a corresponding number of MC simulations, each with ∼10 6 events.We estimate the statistical uncertainty on the A i 's to be at the 10 −2 level.This is less accurate than one could naively deduce due to cancellations that occur when extracting some of the coefficients, but still much smaller than the theoretical error from PDFs and higher-order QCD corrections (see footnote 20 and Ref. [7]).When we consider different m hh categories later in the analysis, the fit is individually performed in each m hh bin through a similar procedure, thus maintaining the same level of statistical uncertainty.
They were computed with our dedicated C þ þ code linked to QCDLOOP [21] and to the LAHPDF routines [22], by setting m h ¼ 125 GeV, m t ¼ 173 GeV and using the CTEQ6LL parton distribution functions.The factorization and renormalization scales have been fixed to Besides m hh , the other kinematic variable that characterizes the gg → hh events is the angle θ between either of the two Higgses and the beam axis in the center-of-mass frame.By total angular momentum conservation, the scattering amplitude can be decomposed into two terms, M 0 and M 2 , describing transitions with, respectively, J z ¼ 0 and J z ¼ AE2, where J z is the projection of J along the beam axis.The amplitude M 2 receives a contribution only from the box diagram and from the dimension-8 operator O gD2 (through the last diagram of Fig. 2).All the other diagrams, instead, mediate J z ¼ 0 transitions.The explicit expression of M 0 and M 2 is reported in Appendix A for convenience.While transitions with all integer values of the total angular momentum J can occur, the leading contributions to M 0 and M 2 come, respectively, from J ¼ 0 (s wave) and J ¼ 2 (d wave).A simple angular momentum decomposition then shows that, up to small corrections, the angular dependence is M 0 ∼ const, M 2 ∼ sin 2 θ. 10 In the SM, the contribution of M 2 to the cross section is always small, as illustrated by Fig. 3: it is negligible at the peak of the m hh distribution (m hh ∼ 400 GeV) and smaller than ∼20% on its tail (m hh ∼ 700 GeV).A shift in the top Yukawa coupling due to NP modifies the value of the box diagram, but cannot change the above conclusion (unless extreme shifts are considered).It is thus interesting to ask whether the dimension-8 operator O gD2 can give a sizable contribution to M 2 through the last diagram of Fig. 2. If this were the case, one could reach a better sensitivity on O gD2 by performing a suitable analysis of the angular distributions of the Higgs decay products.Unfortunately, we find that the effect is numerically small.This is illustrated in Fig. 4, where we show the isocurves of r ¼ σ gD2 =σ tot in the plane ðm hh ; cos θ min Þ.We defined r to be the ratio of the cross section induced by O gD2 alone, σ gD2 , to the total cross section, σ tot , obtained by adding the contribution of O gD2 to the SM.Both cross sections are computed at the partonic level for the process gg → hh and by integrating over angles θ min ≤ θ ≤ π=2.We set the coefficient cgD2 equal to its naïve estimate (8) and choose as benchmark values f ¼ 635 GeV, m Ã ¼ 1.9 TeV, which correspond to g Ã ¼ 3, ξ ≡ ðv=fÞ 2 ¼ 0.15.For r ∼ 0.5 the contribution of O gD2 is as large as the SM one.The plot of Fig. 4 shows that this occurs at the crossover scale m hh ∼ 1.3 TeV, in agreement with the naive estimate of Sec.II B, m hh ∼ y t fðg Ã =y t Þ 1=2 .Increasing the value of θ min (i.e., decreasing cos θ min ) enhances the importance of the J z ¼ 2 component of the cross section, and hence that of O gD2 , but the effect is never large.For example, the value of the crossover scale is reduced at most by ∼10%.We thus conclude that, although an analysis of angular distributions could in principle help disentangle the effects of dimension-8 operators, in practice there is little leverage, and the efficacy of such a strategy is further reduced in the case of the b bγγ final state by the limited range in m hh , which can be realistically probed.In our analysis of gg → hh → b bγγ we will thus make use of the m hh distribution as the main tool to probe the effects induced by NP, neglecting for simplicity angular distributions.

III. DOUBLE HIGGS PHENOMENOLOGY AT 14 TEV AND 100 TEV
The presence of various Higgs decay channels with a non-negligible branching ratio allows the exploration of the Higgs properties in several final states.This is especially true in single Higgs processes, where the relatively large production cross section compensates for the small decay probability in some of the cleanest channels, such as h → γγ and h → ZZ Ã → 4l.In the case of double Higgs production, however, the low signal rate forces one to consider only final states with a sizable branching ratio.This is particularly so at the 14 TeV LHC, where the NLO total SM production rate is around 37 fb −1 , but remains true even at a hypothetical future 100 TeV machine, where the enhanced rates are unavoidably accompanied by larger backgrounds.To obtain a large enough branching ratio, at least one of the two Higgses must decay into a b b pair.For the second Higgs different choices seem possible and have been considered in the literature, namely, h → b b, h → WW Ã , h → τ þ τ − , and h → γγ.The channel hh → b bb b has the highest signal rate (BR ≃ 33.3% in the SM), but is plagued by a large QCD background.Even imposing four b tags, it seems hard to exploit and could require a sophisticated search strategy [3,5,11].The channel with the second highest rate is hh → b bWW Ã , with a branching ratio BR ≃ 24.9% in the SM.Its observation is also threatened by a large background, mainly coming from t t [5,7].It has recently been claimed that a good signal-tobackground ratio can be obtained by using jet substructure techniques and focusing on a very specific region of the parameter space where the two Higgses and their decay products are highly boosted [6].Although encouraging, the results of this analysis suggest that the b bWW Ã final state can be observed at the 14 TeV LHC only with its highluminosity extension.The other channel that has been extensively studied in the literature is hh → b bτ þ τ − [3,5,7,9,13].It is very promising and potentially relevant for the 14 TeV LHC, since its SM branching ratio is sizable, BR ≃ 7.35%, and good signal-to-background ratios seem to be achievable while keeping a relatively large number of signal events.Its actual significance relies, however, on the ability to reconstruct the tau pair, and it will have to be fully assessed by an experimental study.
For the purpose of our study we will focus on the channel hh → γγb b (BR ≃ 0.264% in the SM), which has been considered to be the cleanest one despite its small rate.As shown by previous theoretical studies [4,7,8,10] as well as a recent nonresonant experimental search at ffiffi ffi s p ¼ 8 TeV [23] and a study for the ffiffi ffi s p ¼ 14 TeV high-luminosity LHC [24] by ATLAS, 11 the analysis strategy to observe this decay mode is relatively simple and straightforward.This allows us to avoid unnecessary complications and to concentrate primarily on the interpretation of the search in the context of the EFT description.Clearly, exploring the other available channels from a similar perspective could improve significantly the sensitivity.We leave this interesting follow-up for a future study.Before proceeding with the details of our analysis it is useful to briefly summarize the properties of the kinematic distributions and discuss the main differences between the 14 TeV LHC and a future 100 TeV collider.
It is well known (see, for example, Ref. [12]) that in the SM a cancellation between the box diagram and the triangle diagram involving the Higgs trilinear coupling leads to a depletion of the signal at threshold (m hh ≳ 250 GeV).The peak of the distribution is essentially determined by the fast decrease of the gluon parton distribution functions and is located around m hh ≃ 400 GeV.This conclusion is valid independently of the collider center-of-mass energy ffiffi ffi s p .The main difference between the LHC and a higher-energy collider is a rescaling of the overall cross section (the SM cross section at a 100 TeV machine is around 40 times bigger than the one at the 14 TeV LHC), with small effects on the m hh distribution.As it can be seen from Fig. 5, this latter is modified significantly only on its tail (m hh ≳ 700 GeV), which is enhanced at larger collider energies due to the higher luminosity of gluons.In fact, as we will discuss in Sec.IVA, an important change is also present in the pseudorapidity distribution of the Higgs bosons, which is related to the amount of boost determined by the initial parton energies.
Because of the above considerations, one would naively expect that an analysis similar to the one designed for 14 TeV would continue to work well at 100 TeV.This is, in fact, only partially true.Indeed, the enhanced signal rates can impact the analysis results in two ways.First, they will improve the sensitivity on the parameters, such as the Higgs trilinear coupling, by simply reducing the statistical uncertainty even if the analysis strategy is not modified.Second, 11 See Ref. [25] and Refs.[23,26] for, respectively, theoretical and experimental analyses of the resonant di-Higgs production in the b bγγ final state.because of the larger cross section and enhanced tail, much higher values in m hh will be accessible, which are potentially more sensitive to higher-order operators growing with the energy.For this reason, fully exploiting the potential of a 100 TeV collider requires some modification of our 14 TeV analysis strategy in order to better reconstruct events with a higher boost.Although we will not present a fully optimized analysis for the 100 TeV, we will give a first assessment of how much our final results can improve when jet substructure techniques are used.
As a first step in this direction, it is useful to get an idea of the reach that can be achieved on m hh and p T ðhÞ at different colliders and for different search channels.A complete assessment of this point would require specifying completely the scenario we are interested in and the size of the possible NP contributions to the cross section.To get a rough approximation, however, we can estimate the reach by demanding that a few events (we choose 5 for the estimate) are still present in the tails of the SM distributions above the considered m hh [or p T ðhÞ] value.To be more The cross sections have been computed at NNLO using the k factors of Eq. ( 17).realistic we also assume a 10% overall signal efficiency due to kinematic cuts.The estimates can easily be extracted from the plots of Fig. 6, which show the integrated differential cross sections for the SM signal pp → hh as a function of the lower cut on m hh and of the minimum p T ðhÞ.By including the branching ratio of the most important decay modes we obtain the results shown in Table II.
The reach on m hh is rather limited in the b bγγ channel due to the small signal rate, and even at a 100 TeV collider it does not extend much beyond ∼1.5 TeV.Other channels with larger rates will be crucial to push further the exploration of higher invariant masses.Similar considerations apply also for the reach on p T ðhÞ.The maximal value of p T ðhÞ can be used to estimate whether a given search can benefit from jet substructure techniques or a simple jet analysis is enough.The angular separation between two partons coming from the Higgs decay scales like ΔR ∼ 2m h =p T ðhÞ.To resolve the jets with the standard techniques, we must demand ΔR ≳ R ATLAS;CMS , where R ATLAS ¼ 0.4 and R CMS ¼ 0.5 are the reconstruction cones used by the ATLAS and CMS Collaborations, respectively.It is easy to see that for p T ≳ 500-600 GeV the two partons are not often resolved as two separated objects and jet substructure techniques are likely to be helpful.The results in Table II show that the size of the boosted phase space at 14 TeV is almost negligible for the γγb b channel and very limited for b bb b, γγWW Ã , and b bτ þ τ − .The story changes, however, for an 100 TeV collider.In this case the SM signal can be probed up to m hh ∼ 1.5 TeV in the γγb b channel and up to m hh ∼ 4 TeV for b bb b.In these kinematic regions many boosted Higgses are produced and jet substructure techniques are crucial to reconstruct them.We will analyze this aspect in more detail in Sec.IVA.

IV. ANALYSIS OF pp → hh → γγb b
After discussing the general properties of the pp → gg → hh process, we now focus on the b bγγ final state and describe our analysis strategy.We generated the signal events through our own dedicated code described in Sec.II C. The code includes the 1-loop diagrams of Fig. 2 and is based on the parametrization of Eq. ( 4).It is available from the authors upon request.By using the CTEQ6LL PDFs (LO PDF with LO α s ), setting the factorization and renormalization scales to Q ¼ m hh and m h ¼ 125 GeV, m t ¼ 173 GeV, we find the following LO cross sections for the Standard Model: 16.2 fb and 873.6 fb for 14 TeV and 100 TeV collider energies, respectively.To (partially) include the NLO and NNLO corrections, we rescale the cross sections by the k factors which were computed for the SM in Ref. [27]. 12 The main backgrounds that we considered are the nonresonant processes b bγγ and γγjj and the resonant processes b bh, Zh, and t th (with the subsequent decays h → γγ, Z → b b, and t t → b b þ X).Further backgrounds involving fake photons from jets have been neglected as their estimate is beyond the scope of this work.Experimental analyses of single Higgs production have shown that similar processes in that case can be safely reduced to a subleading level, and we thus assume that this will be possible for double Higgs production as well.We generated all the backgrounds with MADGRAPH5_aMC@NLO v2.1.1 [33] by switching off virtual corrections (i.e., working in LO mode).The output has been interfaced with PYTHIA v6.4 [34] for parton showering and hadronization, and with FASTJET v3.0.6 [35] for jet clustering.The factorization and renormalization scales have been set to the default dynamic scale in MADGRAPH5.Further details about the generation can be found in Appendix B 1. The backgrounds b bγγ, b bh, and Zh have been matched up to one additional jet via the k T -jet MLM matching [36] to partially account for NLO effects.
TABLE II.Lower limits on the pp → hh cross section in the SM and upper limits of the phase space that guarantee 5 signal events.We assume an integrated luminosity L ¼ 3000 fb −1 at the 14 TeV LHC and a 10% efficiency due to kinematic cuts.The numbers in parentheses refer to a 100 TeV collider, with the same integrated luminosity.The cross sections have been computed at NNLO using the k factors of Eq. (17) Compared to the NLO calculation (performed in Ref. [28]), the NNLO corrections give a ∼20% enhancement of the total cross section.Notice that both the NLO and the NNLO corrections to the total rate have been computed in the infinite top mass approximation.It is well known that at LO this approximation leads to very inaccurate kinematic distributions, in particular the m hh one.For this reason the k factors used in our analysis must be considered only as a crude approximation.Recently, a step forward toward a fully differential NLO calculation has been done in Ref. [29], where real emission diagrams were computed at 1-loop, while the finite part of the 2-loop virtual corrections was extracted by resorting to the infinite top mass approximation.For an estimate of the finitetop mass effects and of the dependence of the k factor on a cut on ffiffi ffi ŝ p , see Refs.[30,31].See also Ref. [32] for the resummation of threshold effects in the SM.
In the case of b bγγ we found that allowing for the extra jet increases the cross section by a factor ∼2.We found that this is in good agreement with a complete NLO computation performed with MADGRAPH5_aMC@NLO, indicating that real emissions in this case represent the bulk of the NLO correction.Considering that b bγγ is the dominant background after all cuts, including this effect is very important, but it has been omitted in most of the previous studies.More details about our generation of b bγγ and a thorough discussion of the NLO correction are given in Appendix B 2. For the resonant t th background we apply instead a simple k factor rescaling of the LO cross section to the NLO value reported in Ref. [37].Finally, the γγjj background was simulated without matching and a k factor k γγjj ¼ 2 was applied (this is a more conservative choice than the one proposed in [33] based on an NLO simulation).
To estimate the sensitivity of our analysis we should, ideally, take into account parton-shower/hadronization and smearing in every point of our parameter space.This procedure, however, would be computationally too expensive.Instead we adopt the following simplified approach.We fully extract the signal rate after cuts only for the SM point by performing a hadron-level analysis.For the other points in the parameter space, we apply the same type of analysis to the parton-level samples, and we rescale the signal by the hadron-to-parton cross section ratio computed for the SM, We use the fact that such a ratio is approximately the same both for the SM and in a generic point of the Beyond the Standard Model (BSM) parameter space.The rescaling is performed individually for each of the bins in m hh of Tables V, VII, and VIII shown below.Our analysis correctly takes into account possible nonuniversal signal effects coming from distortions of the kinematic distributions at the partonic level.The rescaling then approximately includes the effect of showering and hadronization.For the jet substructure analysis of Sec.IVA, on the other hand, we follow a different procedure, since there is no parton-level cross section we can use after cuts.We thus compute the BSM signal by starting with the SM one computed at the hadron level and multiplying by the ratio of BSM over SM cross sections before cuts at the parton level This ratio is expected to be approximately the same at the partonic level before cuts and at the hadronic level after cuts if one considers sufficiently narrow bins in m hh , since to very good accuracy the latter is the only variable that controls the kinematics of the signal.We checked that this procedure is accurate at ∼a few% level in the high-m hh categories used for the boosted analysis.
We designed different strategies for 14 TeV and 100 TeV colliders.We begin by describing the one at 14 TeV.Events are triggered by demanding exactly two isolated photons satisfying the minimal reconstruction requirements To forbid extraneous leptonic activity we veto events with isolated leptons (electrons or muons) satisfying Similar isolation criteria are applied to both photons and leptons.A photon (lepton) is considered isolated if the surrounding hadronic activity within a cone of size R ¼ 0.4 (0.3) satisfies p T ðγÞ=ðp T ðγÞ þ p T ðconeÞÞ > 0.9 (0.85).We have checked that more sophisticated photon isolation criteria, like the one proposed in [38], give similar results.
We then cluster the events into R ¼ 0.5 anti-k T jets [39].
We accept only jets with and require that at least two of them are b tagged (the leading two b jets are selected if more than two b jets exist).We assume 70% efficiency for b tagging ϵ b ¼ 0.7, corresponding to 1% of mistag rate ϵ j→b ¼ 0.01 [40][41][42], 13 and 80% efficiency for photon tagging ϵ γ ¼ 0.8 [42,43]. 14After this step, an event consists of two isolated photons, two b-tagged jets, and possible additional jets (whether b tagged or not).
Once the reconstruction of the basic objects by the above procedure is done, we apply the set of cuts described in the following.We first restrict the events to those with two hard photons and two b-tagged jets satisfying  More specifically, we tag with 70% probability only those jets that contain a b hadron with transverse momentum larger than 5 GeV=R jet .For R jet ¼ 0.5 this translates into a minimum 10 GeV transverse momentum.When the jet substructure is used (see Sec. IVA), R jet is set to ΔRðbbÞ, corresponding to the distance between two leading subjets after BDRS declustering.14 Typical rejection rates for fake photons from jets at this working point are in the range 0.1%-0.5%.We will not use this number, however, since backgrounds from fake photons are not included in this analysis.
In the signal, the b b and γγ subsystems tend to be in a back-to-back configuration, and the angular distance between two photons (and similarly between two b jets) is of order ∼2m h =p T ðhÞ ∼ Oð1Þ in the majority of the phase space.While the γγ subsystem in resonant backgrounds has a kinematics similar to the signal, the different origin of the b b pair in each process can be used to distinguish them from the signal.The angular distributions of the signal and of the two dominant backgrounds are shown in Fig. 7.We find that the following cuts (indicated by the dashed lines in Fig. 7), ΔRðb;bÞ < 2; ΔRðγ;γÞ < 2; ΔRðb;γÞ > 1.5; ð24Þ can efficiently reduce the background, especially b bγγ, while retaining most of the signal. 15As a final cut, we restrict the invariant masses of two b jets and of the two photons to the Higgs mass window, 16   105 The resulting cut flow is shown in Table III.We report in Table IV a fit of the signal rate [r ¼ σ × BRðhh → b bγγÞ] obtained after all cuts based on a parametrization analog to that given in Eq. ( 16) for the cross section.At leading order the b bγγ process is initiated by q q and gg and proceeds through diagrams where the two b's tend to emerge with a large relative angle.Diagrams initiated by gq and g q at next-to-leading order lead to topologies that can more easily fake the angular configuration of the signal but account for only a minor fraction of events.16 The width of the interval 120 < m reco γγ < 130 GeV corresponds to ∼3 times the experimental resolution on photon pairs; see [44,45].The same mass window was adopted by the recent CMS study of di-Higgs resonant production [26]; see also Ref. [23].The width of the interval 105 < m reco b b < 145 GeV corresponds to ∼2 times the experimental resolution on b-jet pairs from Higgs decays (after correcting for resolution effects); see [46,47].We do not smear photons in our simulation, nor do we include a specific efficiency for the reconstruction of the photon pair.The mass window on m reco γγ is, on the other hand, sufficiently wide that almost all of the signal is retained, so that the efficiency would be close to 100% even including a finite energy resolution on photons.The efficiency of the cuts of Eq. ( 25) reported in Table III corresponds, in the case of the signal, to the efficiency for the reconstruction of the b b pair (equal to ∼0.72).
With our analysis strategy, γγb b and t th are the two major backgrounds.The latter tends to produce extra jets from the hadronic decay of the W's.One could thus consider applying a veto on the extra hadronic activity to enhance the signal significance.The potential impact of a jet veto can be seen in Fig. 8.For example, further restricting the events to the region with NðjetsÞ < 4, in addition to all the previous cuts, can remove roughly 80% of the t th background while keeping ∼70% of the signal.Alternatively, one could look for hadronic W's by iteratively forming jet pairs with invariant mass lying in a given window around the W mass; for example, we will use (70,100) GeV in the following.As illustrated in Fig. 8, vetoing hadronic W's in addition to our cuts can reduce ∼50% of the t th background while retaining more than 90% of the signal.We find that applying either of these cuts leads to a modest increase in the significance, at the cost of reducing the number of signal events.We thus decided not to exploit any form of extra-jet vetoing at 14 TeV, motivated by the necessity of retaining as many signal events as possible.
As discussed in Sec.II C, several diagrams with different energy scalings contribute to the signal.We can thus use the invariant mass of the reconstructed hh system to differentiate the various effects and improve the sensitivity on the Higgs couplings.To this purpose the events are subdivided into six different categories in m reco hh .The corresponding numbers are reported in Table V and will be used in Sec.V to extract our bounds on the coefficients of the effective operators.
Let us now discuss the case of a 100 TeV collider.We start by considering a strategy similar to the one adopted at 14 TeV and different from this one mostly for the numerical values of the cuts and selection parameters.The nominal p T threshold on the reconstructed jets is increased to p T ðjÞ > 35 GeV to take into account the busier environment of the higher energy collision.We define the same pseudorapidity region for the acceptance of jets, photons, and leptons although more signal events fall into the high-jηj region due to the boost along the beam axes (this point will be discussed in detail in Sec.IVA).Similarly, in our initial Jet multiplicity After the above cuts we find angular distributions similar to those shown in Fig. 7, and we thus apply the same cuts as in Eqs.(24) and (25).The cut flow at ffiffi ffi s p ¼ 100 TeV is reported in Table VI.Different from the 14 TeV case, at this level the dominant background is t th, rather than b bγγ, due to its steeper increase with the collider energy.Imposing a veto on extra hadronic activity is beneficial to reduce this background without strongly affecting the signal, as illustrated by the distributions of Fig. 9.We find that a veto on hadronic W's is most efficient to maximize the signal significance, and thus apply it.The final number of signal and background events after the W veto is reported in the last row of Table VI, while a fit to the total signal rate after all cuts is given in Table IV.As a final step, we subdivide events into six m reco hh categories and report the corresponding numbers in Table VII.

A. Recovering the boosted topologies
As we anticipated in the previous discussion, a possible issue with our analysis strategy is the loss of sensitivity for the kinematical configurations containing boosted Higgses.Although this effect is not likely to be relevant for the 14 TeV LHC, it can potentially affect the gg → hh searches at a future higher-energy collider.In this section we present a first estimate of the improvement that can be achieved in the reconstruction of highly energetic events by the use of jet substructure techniques.Notice that the development of a fully optimized analysis will most likely require a hybrid strategy that smoothly interpolates between a traditional jet analysis and one using boosted techniques [48].Devising such a strategy is beyond the scope of the present work, and we postpone it to a future study.In this section we instead adopt a "minimal" approach and use either the standard analysis or the jet substructure technique in each m reco hh bin.There are basically two different effects that lead to boosted Higgses, namely, the boost of the whole hh system along the beam axis and the production of a Higgs with high p T .Both effects will become relevant at a future highenergy collider.At 100 TeV the di-Higgs system acquires on  average a non-negligible amount of momentum along the beam axis, and, as a result, the events more frequently leak into the high-jηj region.The situation is illustrated by the left plot of Fig. 10, which suggests that a larger pseudorapidity range for object reconstruction and flavor tagging could improve the sensitivity to the gg → hh process.For instance, at the 14 TeV LHC roughly 13% of the partonic signal lays outside the acceptance region jηj < 2.5 (solid vertical line in the left panel of Fig. 10).If one assumes the same acceptance region, this fraction increases to ∼30% at a 100 TeV collider.To obtain the same acceptance efficiency of the 14 TeV LHC, the coverage region should be extended to jηj < 3.3 at the higher-energy collider (vertical dotdashed line in the left panel of Fig. 10). 17 A further problem with boosted events is that our analysis relying on a dijet-style search fails when the two b jets are not well separated.A similar argument can be applied to the diphoton system.The boosted Higgses are more likely to be produced in the high invariant mass tail of the m hh distribution.This effect is illustrated in the right panel of Fig. 10, where we plot the ΔR separation between two b's or photons at the partonic level for the two categories of events with m hh > 1000 GeV and m hh < 1000 GeV.Most events belonging to the lower m hh category have a clear two-prong topology, whereas a significant number of events in the higher m hh category fail to be resolved as two well separated partons.
For our present purposes it will be sufficient to perform a simple analysis on the b b subsystem and compare the performance of our traditional jet-based analysis with the substructure method.For the γγ subsystem we consider two possibilities: a reduced isolation cone size (R iso ¼ 0.2) and an ad hoc mutual photon isolation (a somewhat similar prescription was proposed in [49] for a boosted di-tau system).In the latter, photon isolation is imposed neglecting the other photon in the γγ system, and an event is retained as long as ΔRðγ; γÞ > 0.2.This procedure is quite similar to our traditional jet-based analysis up to the slight modification of the photon isolation criterion.We do not modify instead the isolation criterion on leptons. 18 After this step, the events are clustered into R ¼ 1.5 "fat jets" with the Cambridge/Aachen jet algorithm [50,51].We iteratively look over those fat jets and apply the BDRS subjet-finding technique of Ref. [52] (with the same declustering parameters as in [52]).To ensure enough boost we only look into fat jets with p T ðjÞ > 150 GeV.The declustering algorithm stops when it successfully identifies three subjets (the third subjet takes into account the leading gluon emission from either bottom quark line).The fat jet is identified as a Higgs jet if the invariant mass of the three subjets falls into the Higgs mass window m reco b b ¼ ½105; 145 GeV [see Eq. ( 25)].If multiple candidates exist, we pick the one whose invariant mass is closest to the Higgs mass.Two b taggings are performed at subjet level assuming a 70% b-tagging efficiency.The cuts on the angular separations in Eq. ( 24) are applied to the subjets and the photons.Finally, we do not apply any veto on extra hadronic activity or on hadronic W bosons.
The performances of the substructure analysis on signal events are summarized in Table VIII and compared to the traditional jet-based analysis.We consider three possible scenarios that differ by the treatment of the γγ system: "Substructure I" uses a standard photon isolation criterion with R iso ¼ 0.4, whereas in "Substructure II" the cone size is reduced to R iso ¼ 0.2."Substructure III" corresponds instead to the mutual isolation criterion.Not surprisingly, the substructure analysis becomes more efficient compared to the traditional one for m reco hh ≳ 1000 GeV; this matches our naive expectation based on the ΔRðb; bÞ and ΔRðγ; γÞ distributions of Fig. 10.For m reco hh ≲ 1000 GeV, on the other hand, only part of the events is correctly reconstructed using the substructure technique (due to the finite fat-jet size R ¼ 1.5), and the traditional analysis is more efficient.This is also clearly illustrated by the plots of Fig. 11, which show the m reco hh distributions for signal events with both After imposing the cuts in Eqs. ( 23) and ( 24) on partonic signal events, these fractions reduce to ∼6.5% for the 14 TeV LHC and ∼24% for the 100 TeV case.To recover the 14 TeV acceptance efficiency at the higher-energy collider, the pseudorapidity region should be extended to jηj < 3.6.

18
While naively one would expect that leptons from the boosted tops in the resonant t th background fail more frequently the isolation requirement, we find that the efficiency drop is not significant for the m hh range of interest.
strategies.These results suggest that the use of jet substructure in the b bγγ final state is crucial only at a 100 TeV collider, while it has a minor impact at the 14 TeV LHC, where large values of m hh are suppressed by the fast drop of the gluon luminosity.The relevance of jet substructure can be determined in a more quantitative way by comparing the reach in m reco hh at 14 TeV and 100 TeV.At the partonic level, Fig. 6 shows that the same number of events one has at 14 TeV after imposing the cut m hh > 850 GeV is obtained at 100 TeV for m hh > 2570 GeV.At the hadronic level, the increase in the reach is much more modest than the naive expectation when adopting a traditional jet-based analysis.For instance, the signal rate at 14 TeV is N events ∼ 0.5 for m reco hh > 850 GeV.The same number of events is obtained at 100 TeV for m reco hh > 1540 GeV using a traditional analysis, and m reco hh > 2560 GeV using jet substructure with the reduced photon isolation cone R iso ¼ 0.2.Hence, a jet substructure analysis (along with the appropriate modification of the photon isolation criteria) is essential to get close to the naive expectation.
Focusing on the 100 TeV case, we show in Table IX the number of signal and background events expected with the  substructure analysis and a reduced photon isolation cone size (scenario Substructure II) in six categories with high m reco hh .We include only the dominant b bγγ and t th backgrounds for simplicity. 19A further reduction of the t th background could in principle be possible by iteratively looking for boosted top and W jets through dedicated tagging algorithms, and by discarding events where such objects are reconstructed.This is, however, beyond the scope of the present work.On the other hand, we have checked that a simple veto on hadronic W's reconstructed from pairs of jets, as applied in the traditional analysis, is not efficient on events with such high m reco hh .The results in Table IX will be used in Sec.V to assess the impact of jet substructure in setting bounds on the coefficients of the effective operators.

V. SENSITIVITY ON THE EFT COEFFICIENTS
The results of the previous section have been processed through a Bayesian statistical analysis to extract the sensitivity on the coefficients of the Higgs effective Lagrangian.In each case, a probability distribution of the relevant parameters is obtained by constructing a likelihood function from the signal and background number of events in the six m hh categories, and marginalizing over (or fixing) the remaining parameters.The injected signal is the SM one for all the results presented in this section.A flat prior is assumed for the coefficients of the Higgs effective Lagrangian, except when they are constrained by single-Higgs measurements.In the latter case, we use the ATLAS projections for the highluminosity LHC of Ref. [53] to construct a likelihood function and use it as a prior.Reference [53] reports the estimated precision on various Higgs decay channels expected for ffiffi ffi s p ¼ 14 TeV with integrated luminosities L ¼ 300 fb −1 and L ¼ 3 ab −1 .We will further study the case of a future proton-proton circular collider operating at ffiffi ffi s p ¼ 100 TeV with 3 ab −1 of integrated luminosity.When single-Higgs measurements are required to perform marginalization in this latter scenario, we use the ATLAS projections for L ¼ 3 ab −1 at ffiffi ffi s p ¼ 14 TeV, lacking a specific estimate of the precision reachable at a 100 TeV machine.We will thus consider the following three benchmark scenarios: For simplicity, we do not include in our analysis the theoretical uncertainty on the prediction of the signal cross section nor the systematic uncertainties that will characterize the extraction of the background from data in a realistic experimental analysis. 20GeV) The solid line is obtained by adopting the substructure technique with reduced photon isolation cone size (scenario Substructure II described in the text), whereas the dot-dashed line corresponds to a traditional jet-based analysis.19 We have checked that Zh and b bh are negligible, as the corresponding number of events in each category is always smaller than 0.1.We expect that also the background γγjj can be reduced to a subdominant level as long as the efficiency for making two b tags at the subjet level is sufficiently high.The γγb b background was newly generated for the substructure analysis with relaxed generation cuts ΔRðb; bÞ > 0.1, ΔRðγ; γÞ > 0.15; these are less stringent than those specified in Appendix B 2.
As a first result, we derive the precision that can be obtained on the signal strength multiplier μ ¼ σ=σ SM in the three benchmark scenarios described above.No marginalization is performed in this simple case.We obtain the following 68% intervals on μ: Very similar results are obtained from an inclusive analysis without m hh categories, as one naively expects since μ accounts only for an overall rescaling of the total cross section of the signal.
Turning to the determination of the coefficients of the effective Lagrangian, we first consider the case of the nonlinear parametrization of Eq. ( 4). Figure 12 shows the 68% probability contours in the planes ðc 2t ; c 3 Þ and ðc 2t ; c 2g Þ obtained through the procedure described above.In this case the marginalization is performed over two parameters: c t , with a prior obtained from single-Higgs measurements, and the branching ratio for the decay hh → b bγγ.For the latter parameter we assume a Gaussian distribution around the SM value with standard deviation equal to 0.15 and 0.10, respectively, for L ¼ 300 fb −1 and L ¼ 3 ab −1 (at both 14 TeV and 100 TeV).For simplicity, the remaining couplings are set to their SM values: c g ¼ c 2g ¼ 0 in the plot on the left of Fig. 12; c 3 ¼ 1 and c g ¼ 0 in the plot on the right.We have checked that performing an additional marginalization over c g slightly decreases the precision on the measured couplings, without changing the shape of the contours.
We find that the couplings c 3 and c 2t are strongly anticorrelated, and the precision expected on c 2t is much higher than the one on the Higgs trilinear coupling c 3 .This is in agreement with previous studies, which pointed out the strong sensitivity of the double Higgs cross section on the t thh quartic coupling; see [12] and references therein.The coupling c 2g is determined even more accurately, although its naive estimate is suppressed, compared to that of c 3 and c 2t , by two powers of a weak spurion if the Higgs is a pseudo-NG boson; see Eqs. ( 5) and ( 3).Although the Higgs trilinear coupling c 3 is the less accurately measured parameter, its precision can be highly increased at a 100 TeV collider, as shown by the left plot of Fig. 12.This is mainly due to the higher number of signal events that can be produced at this machine.Using boosted jet techniques does not seem to improve significantly the accuracy on c 3 and c 2t , while it increases dramatically that on c 2g .This is expected, since jet substructure techniques  are most relevant to reconstruct highly boosted events with large m hh , which are especially important to determine c 2g .The situation is illustrated by Fig. 13, where the red dotted curve is obtained by including the first five categories (250 GeV < m hh < 1000 GeV) of the traditional jet-based analysis and the last five categories (m hh > 1000 GeV) of the Substructure II analysis; see Tables VII and IX.
To break the degeneracy among the various parameters and extract them precisely, it is crucial to make use of the information on m hh .We find that an inclusive analysis, where events are not classified in the six m hh categories of Tables V and VII, is much less powerful in constraining the Higgs couplings, especially in the case of a 100 TeV collider.This is illustrated by the plots of Fig. 14 in the plane ðc 2t ; c 3 Þ.It is evident how at 100 TeV only an exclusive analysis can break the degeneracy between c 2t and c 3 .As expected from the discussion of Section II C, categories with larger m hh most strongly constrain c 2t and c 2g , while c 3 is mainly determined by events at threshold; see Fig. 15.Given the relevance of the categories with large m hh in determining the Higgs couplings, it is important to give an assessment on the validity of the EFT approach in our analysis.Following the discussion of Sec.II B, we evaluate the minimum coupling strength g min ∼ ffiffi ffi δ p ð Ē=vÞ by estimating the maximal invariant mass Ē of the events that lead to a precision δ on c 2g (i.e., the parameter controlling the term that grows quadratically with the energy in the scattering amplitude).To this aim we rederive the 68% probability contours in the plane ðc 2t ; c 2g Þ by removing one or more of the categories with large m hh , as a way to estimate their impact on the determination of c 2g .We show such plots in Fig. 16 for the HL-LHC and the FCC 100 .We find LHC 14 HL-LHC FCC 100 FCC 100 (boosted) where the last column refers to the analysis including jet substructure at the FCC 100 (first five categories of traditional analysis plus last five categories of Substructure II analysis).As discussed in Sec.II B, the value of g min determines the extension of the region that can be probed under the validity of the EFT; see Fig. 1.In particular, g min < y t only in the case of a 100 TeV collider, which suggests that an analysis including only dimension-6 operators at the LHC (even at its high-luminosity upgrade) is not sensitive to the case of a pseudo-NGB Higgs with fully or partially composite t R , although it can probe the case of a generic composite Higgs.FIG.16 (color online).The 68% probability contours in the plane ðc 2t ; c 2g Þ for the HL-LHC (left plot) and the FCC 100 (right plot).The different curves have been obtained by removing the following m hh categories of Tables V and VII from the fit: 6 (dashed red line); 6 and 5 (dot-dashed brown line); 6, 5, and 4 (dotted blue line).The continuous black contour is obtained by including all the categories in the fit.
We now turn to the effective Lagrangian for a Higgs doublet, Eq. ( 2), and show the constraints on its coefficients.The production cross section, in this case, depends on cH , cu , cg , and c6 , while the branching ratio for hh → b bγγ depends on cH , cu , cg as well as on cd and cγ .The latter two coefficients parametrize, respectively, the modification of the down quark Yukawa couplings, in particular that of the bottom quark, and the contact interaction between the Higgs boson and two photons.Their precise definition can be found in Ref. [15], and it is analogous to that of cu and cg in Eq. (2).We will set cγ ¼ 0 for simplicity in the following.
We start considering the constraints on cu and cg that follow from both single and double Higgs production.The 68% probability contours are shown in Fig. 17 for the three benchmark scenarios considered in this section.The solid blue curve refers to double Higgs production in the b bγγ final state (our analysis), while the dashed and dotted blue curves correspond to the constraint from, respectively, all single-Higgs processes except t th, and t th alone 21 (the ATLAS projections from Ref. [53]).They have been obtained by fixing cH ¼ cd ¼ 0 and marginalizing over c6 with flat prior (in the case of double Higgs production).The two dashed ellipses in the plots in the upper row correspond to the two degenerate solutions that follow from the interference between cg and cu in the gg → h cross section.The filled areas denote instead the 68% probability regions obtained by combining the following processes: all single Higgs processes including t th (orange region); all FIG.17 (color online).The 68% probability contours in the plane ðc u ; cg Þ for the three benchmark scenarios: LHC 14 (upper left plot); HL-LHC (upper right and lower left plots); FCC 100 (lower right plot).Different curves correspond to double Higgs production in the b bγγ final state (solid blue line); all single-Higgs processes except t th (dashed blue line); t th alone (dotted blue line).Filled areas correspond to all single Higgs processes including t th (orange region); all single Higgs processes plus double Higgs production (dark blue region); all single Higgs processes except t th plus double Higgs production (light blue area in the lower left plot).

21
Notice that the dotted curves are not vertical in the plane ðc u ; cg Þ due to the dependence of the total Higgs width on cg .A further dependence on cg would follow from the contribution of diagrams with an insertion of O g to the t th cross section; this effect has not been included in our fit.We thank C. Grojean for drawing our attention to this point.
single Higgs processes plus double Higgs production (dark blue region); all single Higgs processes except t th plus double Higgs production (light blue area in the bottom left plot of Fig. 17).They have been obtained by marginalizing over cH , cd , and c6 .One can see that double Higgs production breaks the degeneracy on cg that remains even after including t th among single-Higgs processes.It also helps increase the precision on both cu and cg compared to that following from single Higgs processes alone, especially at a 100 TeV collider.In this respect, it is interesting to notice that cu and cg are strongly correlated when considering single-Higgs measurements without t th, and that the extension of the dashed ellipses is much larger than the error obtained on each individual parameter by setting the other to its SM value (in fact, the ellipses become infinite bands if one lets cγ vary, as a consequence of the degeneracy in the h → γγ decay rate). 22The process t th mainly constrains cu and is crucial to reduce the overall experimental uncertainty.In comparison, gg → hh → b bγγ is less powerful but still competitive in constraining cu .For example, by removing t th from the combination of single and double Higgs at the HL-LHC, one obtains the light blue region (instead of the orange one) that is sensibly smaller than the dashed contour.One must also notice that the t th dotted curve in Fig. 17 corresponds to the combination of several decay channels (all those studied by Ref. [53]: h → γγ; ZZ; μμ), while the continuous blue curve denotes double Higgs production in the b bγγ final state alone, i.e., the one studied in this work.Combining with other final states (e.g., b bττ and b bWW) will make double Higgs production even more competitive for the determination of cu , as well as of cg .
The plot on the left of Fig. 18 shows the 68% probability contours in the plane ðc u ; c6 Þ, obtained by marginalizing over cH , cd , and cg .As already suggested by Fig. 12, the precision on c6 (i.e., the parameter that controls the modification of the Higgs trilinear coupling) is much smaller than that on cu .At ffiffi ffi s p ¼ 14 TeV, in particular, the constraint is on values of c6 larger than 1. 23 Further marginalizing over cu leads to the probability functions shown in the right plot of Fig. 18 for the three benchmark scenarios.Notice that even at the HL-LHC the likelihood is far from being Gaussian, and a second maximum is present at c6 ≃ 4.5.We find the following 68% probability intervals on c6 :  By fixing all the parameters to their SM value but one, the likelihood obtained from Ref. [53] is approximately Gaussian.By including all single Higgs processes except t th we obtain the following standard deviations: 23 See the discussion of Sec.II A on the validity of the effective field description in this case. of Fig. 19 (dotted red curve).The two solutions are in this case almost degenerate, and, in fact, the second (non-SM) maximum is slightly higher than the one at c6 ¼ 0 as an effect of the marginalization on cu . 24A fully exclusive analysis thus removes the degeneracy and helps reduce the significance of the unphysical solution.
Another interesting question concerns the relevance of the marginalization over cH , cd , cg , and cu in determining the Higgs trilinear coupling.The right plot of Fig. 19 shows the probability function that is obtained by switching off the marginalization in the HL-LHC scenario (dotted red curve).The effect of marginalization is that of flattening the SM solution, although in a marginal way at 14 TeV.The effect is, on the other hand, much more important at 100 TeV.On a more quantitative side, the 68% probability intervals that follow without marginalization are ½−1.4;5.8 (LHC 14 ), ½−1.0; 1.6∪½3.8;5.1 (HL-LHC), and ½−0.18; 0.18 (FCC 100 ).It is clear that a precise measurement of the Higgs trilinear coupling, as it follows at 100 TeV from the b bγγ channel or as it might be obtained at the high-luminosity LHC from a combination of several decay modes, relies on the accurate extraction of the other couplings.In particular, a significant uncertainty on c6 could follow from a poor determination of cu .This is illustrated in Fig. 20 for the 100 TeV scenario.The plot shows the extension of the 68% probability interval on c6 as a function of σðc u Þ, where cu is marginalized with a Gaussian distribution with standard deviation σðc u Þ.The uncertainty on c6 grows with σðc u Þ for small values of this latter parameter, while it becomes constant for σðc u Þ ≳ 0.2.In this limit double Higgs production alone determines cu with a smaller error, independently of single-Higgs measurements.
The impact of the uncertainty on the top Yukawa coupling in a precise determination of the Higgs trilinear was already discussed in Ref. [54], although in a scenario where the only modifications with respect to the SM are in the values of these two couplings.Such a scenario of NP, however, is a disfavored one.If the Higgs boson is part of a weak doublet and the new states are heavy, one gets the effective Lagrangian (2) at low energy, where the same operator O u that gives a modification to the top Yukawa coupling also generates the quartic t thh interaction, whose contribution must thus be included.On the other hand, it is This "unwanted" feature can be avoided by profiling, instead of marginalizing, over cu .The second maximum in this case would be lower than the one at c6 ¼ 0, since the highest peak of the two-dimensional likelihood is indeed the one at the SM point.
possible to get a modification of the top Yukawa without having a direct t thh interaction if the new states are light, as, for example, in two-Higgs doublet models.In that case, however, the new states will give an extra contribution to the production cross section, which is, in fact, required to match that from the t thh vertex in the decoupling limit.One can envisage a situation where the top Yukawa coupling c t is modified by new heavy physics while c 2t ¼ 0 if the Higgs does not belong to a weak doublet, as, for example, in the case of a dilaton.This class of models is, however, disfavored by current data.These considerations suggest that any NP that modifies the top Yukawa coupling is likely to have an additional impact on the double Higgs cross section.More specifically, if the new states are heavy and the Higgs is part of a doublet, then any uncertainty on the value of the top Yukawa coupling will reflect on a contribution from the t thh vertex that must be properly included.Of course, one can restrict oneself to the SM case and ignore the t thh anomalous coupling; in that case, however, the top Yukawa coupling can be more precisely determined from the top quark mass, and its uncertainty is sufficiently small to have a negligible impact on double Higgs production.

VI. COMPARISON WITH PREVIOUS WORKS
The results of the previous section suggest that an analysis of double Higgs production in the b bγγ final state can determine the Higgs trilinear coupling with a ∼30% precision at the FCC 100 . 25At the HL-LHC, on the other hand, the estimated uncertainty is at the level of 200%, and, in fact, another solution of the fit is present for a trilinear coupling approximately equal to 5.5 times its SM value.These results are more pessimistic than other estimates that appeared in the literature, and it is thus useful to compare our analysis with previous works in order to highlight the different strategies and assumptions.
The first detailed analysis of the b bγγ final state was performed by Baur, Plehn, and Rainwater in Ref. [4].They retained events with one or more b tags at the LHC, while two b tags were required at its high-luminosity phase.An analysis with only one b tag seems much more involved due to the larger background and the systematic error from combinatorics.We will thus compare the results for the HL-LHC obtained with two b tags.A first difference with our analysis is in the value of the b-tagging efficiency and the corresponding light jet mistag rate, since Ref. [4] assumes ϵ b ¼ 0.5 and ϵ j→b ¼ 1=23, to be compared with ϵ b ¼ 0.7 and ϵ j→b ¼ 0.01 used in the present work (for both the LHC and its high-luminosity upgrade).While the choice of Ref. [4] was based on early MC simulations, more recent results based on LHC data have shown that larger b-tagging efficiencies can be achieved while maintaining acceptable levels of rejection rates [40][41][42].As a consequence of the larger value of ϵ j→b , the background jjγγ is found to be the largest one in Ref. [4], different from our estimates.A second difference concerns the estimate of the background b bγγ, which we found to be enhanced at NLO by a large k factor k b bγγ ≃ 2. Reference [4] instead computes all the backgrounds at LO and rescales their cross sections by an ad hoc common factor 1.3 to account for NLO effects, thus underestimating b bγγ.Finally, a narrower mass window for the photon pair, m reco γγ ¼ m h AE 2.3 GeV, is assumed in the study of Ref. [4].This leads to a reduction by a factor ∼2 compared to our estimates of all the nonresonant backgrounds, since their m reco γγ distribution is nearly flat around the Higgs mass.We have chosen instead to make a more conservative cut on m reco γγ and to leave the issue of optimizing the width of the mass window to a fully realistic experimental analysis.As a consequence of the above different assumptions, the total background in Ref. [4] is rather smaller than in our analysis, which led the authors to a more optimistic estimate of the precision on the trilinear coupling.
A more recent analysis of the b bγγ channel has been carried out in Refs.[7,8,10].Reference [7] by Baglio et al. estimates 47 SM signal events after all cuts with L ¼ 3 ab −1 at the high-luminosity LHC, which should be compared with our ∼13 events.Although our smaller rate can in part be explained by the inclusion of the photon efficiency factor 0.8 2 ¼ 0.64 (not included in Ref. [7]), we were not able to fully identify the origin of such a difference.Furthermore, the b bγγ continuum is found to be negligible in Ref. [7] due to a limited MC statistics, thus leading to an underestimation of the background rate.The comparison with Ref. [10] is instead difficult, since the number of signal and background events is not reported.The same kinematic cuts of Ref. [7] are, however, applied, and it is said that an agreement is found with the efficiencies there reported.This suggests that the same underestimation of the background also plagues Ref. [10].
Finally, we compare with the analysis by Yao in Ref. [8], upon which the Snowmass study [55] is based.This is a sophisticated study including a full-fledged simulation of showering, hadronization, and detector effects.The estimated numbers of SM signal and background events at the LHC after all cuts are, respectively, 16.6 and 53.4 (40.1 of which from b bγγ) for L ¼ 3 ab −1 , with a corresponding statistical significance S= ffiffiffi ffi B p ¼ 2.3.This is compatible with the numbers reported in Table III, which give S= ffiffiffi ffi B p ¼ 2.1.Our analysis thus agrees with the signal and background rates estimated by Yao.Based on these numbers, this latter calculates a precision on the trilinear coupling equal to 50% and 8%, respectively, for the HL-LHC and the FCC 100 .This is much more optimistic than our corresponding estimates 200% and 30% reported 25 Here and in the following, by sensitivity/precision on some Higgs coupling or coefficient of the effective Lagrangian we mean the 68% error on its measured value for an injected SM signal.
in Sec.V.The discrepancy follows because the uncertainty on the trilinear is derived in Ref. [8] from the dependence of the signal cross section on this parameter before cuts.In particular, a linear approximation is used around the SM point with a slope dðσ=σ SM Þ=dc 3 ≃ −0.8 taken from Ref. [7].However, the slope decreases substantially (in absolute value) after imposing all cuts: at the end of our analysis we find dðσ=σ SM Þ=dc 3 ≃ −0.58 and −0.50, respectively, for the HL-LHC and the FCC 100 . 26By using these numbers, the uncertainties on the trilinear coupling that follows from the rates of Ref. [8] are equal to 90% and 16%, respectively, for the HL-LHC and the FCC 100 . 27The difference between these improved values and our results is finally due to the fact that the linear approximation is not accurate for ffiffi ffi s p ¼ 14 TeV (a second solution exists), and to the marginalization on cu , cH , cd , cg (not included in Ref. [8]), which has a large effect at ffiffi ffi s p ¼ 100 TeV as explained in Sec.V.
Besides the above theoretical studies, the ATLAS Collaboration has recently completed a Monte Carlo analysis of the b bγγ decay mode for the HL-LHC [24].The simulation of the b bγγ background is done by matching up to one extra jet at the matrix-element level, thus properly including the bulk of the NLO corrections.The reported signal and background rates are consistent with those obtained in our analysis, although backgrounds with fake photons (b bγj) and fake b jets from charm quarks (ccγγ) are found to be sizable, different from what is assumed in this work.It is interesting to see if further experimental analyses will confirm this preliminary study and strategies will be found to reduce the size of these reducible backgrounds.

VII. SUMMARY AND OUTLOOK
In this work we presented an analysis of double Higgs production via gluon fusion.The novelty of our approach with respect to most of the previous studies is the use of an EFT perspective that allows one to encode all possible effects from heavy NP into a small set of deformations of the SM Lagrangian.Although the leading contributions to observables are in general expected to arise from dimension-6 operators, when these latter are suppressed due to selection rules higher-dimensional operators can become important.We pointed out that this occurs in double Higgs production, where the dimension-6 operator O g violates a Higgs shift symmetry and is thus suppressed if the Higgs is a pseudo-NG boson.In this case dimension-8 operators become relevant in the high invariant-mass tail of the kinematic distributions.A careful assessment of the range of validity of the EFT description is thus required, like the one provided in Sec.II.There we analyzed in detail the implications of the standard SILH power counting, which predicts Oðv 2 =f 2 Þ deviations in a large set of observables, including the Higgs trilinear coupling.We discussed alternative scenarios, including a Higgs-portal model, that imply a modified power counting and lead to large deviations in the trilinear coupling while predicting small modifications to the other Higgs couplings.These presumably will be the first scenarios to be probed by early measurements of double Higgs production at the LHC, where the precision on the trilinear coupling is expected to be limited.
We provided an explicit application of the EFT approach by assessing the experimental sensitivity on the coefficients of the relevant local operators.For definiteness we focused on the gg → hh → γγb b process, which has been recognized as one of the cleanest channels to exploit double Higgs production.Although this process was previously analyzed in several works, our study provides some important improvements that led to significantly different results (see Sec. VI for a full comparison).One of the most relevant observables that can be extracted is the Higgs trilinear coupling c 3 .We found that the γγb b channel allows the determination of c 3 with a fair precision (∼30% with L ¼ 3 ab −1 ) only at a future 100 TeV hadron collider.The prospects for the LHC, instead, are much less optimistic, and only an Oð1Þ determination seems possible even with 3 ab −1 of integrated luminosity.This result is significantly worse than what is usually claimed in the literature.The origin of the discrepancy is in large part due to a more careful determination of the background processes.In particular, we found that the irreducible b bγγ background is enhanced by a sizable NLO k factor (k b bγγ ∼ 2), not included in most of the previous works, and is thus larger than previously thought.We also showed that the determination of c 3 is affected by the choice of the statistical treatment.In particular, a full fit including possible corrections to all relevant dimension-6 operators leads to significant differences with respect to the procedure, often used in the literature, where only the Higgs trilinear coupling is allowed to vary.An important role is played by the uncertainty on the coupling of the Higgs boson to the top quark.The reason is that, in the context of the effective Lagrangian for a Higgs doublet, the same operator O u that controls the modification to the top Yukawa coupling also generates a quartic t thh anomalous interaction that leads to a new contribution to the double Higgs production cross section.This latter effect must be properly included when estimating the precision on the Higgs trilinear coupling.We find that an uncertainty of order 10% on the measurement of cu can nearly double the error on c 3 at a future 100 TeV collider (see Fig. 20).One 26In fact, the dependence on the trilinear coupling also depends on the m hh category considered, the categories with lowest m hh having the largest absolute slope, as expected.

27
During the completion of this paper, an updated analysis by Yao was presented in a conference talk [56], where a similar precision is obtained at 100 TeV after taking into account the correct slope.can, in fact, turn the argument around and use double Higgs production to determine cu .Our results show that the accuracy one can obtain in this way is potentially competitive with the one from the t th process (see Fig. 17), especially if several final states are combined.On more general terms, other couplings can be extracted from the gg → hh process.In particular, producing a pair of Higgses gives unique access to the t thh and gghh quartic interactions, parametrized, respectively, by c 2t and c 2g .These should be regarded as independent parameters from the wider perspective of a nonlinear effective Lagrangian.We found that these couplings can be determined much more accurately than c 3 , as shown in Fig. 12, in agreement with previous works.
Another important point of our study is the effectiveness of an exclusive analysis exploiting the m hh differential distribution.This procedure is particularly useful to disentangle the contributions arising from different effective operators since each of them leads to very specific deformations of the m hh distribution.This analysis strategy can be relevant for the LHC mainly in its highluminosity phase and can have a dramatic impact at a future 100 TeV collider (see Figs. 14 and 15).The other kinematic variable characterizing double Higgs events, namely the angular separation between the two Higgses, was instead found to have a marginal role in our analysis.The expected distribution, indeed, is almost flat in the SM and nearly unchanged by NP effects.A possible exception to this result is the deformations due to dimension-8 operators.These are, however, unlikely to be accessible at the LHC, especially in the rare b bγγ decay mode, although they are expected to be relevant at future colliders.We also explored the possibility of using jet substructure techniques to improve our analysis strategy.We found that these methods can be relevant in the high invariant-mass tail of the m hh distribution, namely, for m hh ≳ 1 TeV.The impact of such improved analysis seems to be marginal at the LHC, whereas at future higher-energy colliders it can lead to a significant increase in the sensitivity on O g (see Fig. 13) and, possibly, on dimension-8 operators.
There are a few directions in which our analysis could be improved, which we leave for future investigation.In the present work, for example, we did not try to optimize the mass windows for the reconstruction of the b b and γγ pairs.The choice of these windows has a strong impact on the nonresonant backgrounds, as they scale roughly linearly with the window size, and it is not unreasonable to expect that some improvement could be achieved with a more sophisticated analysis.Another point that we did not fully explore is the estimate of the backgrounds with fake photons.The analogy with single-Higgs production suggests that it should be possible to reduce these backgrounds to a subdominant level, but a thorough quantitative analysis is required to clarify this issue.
As previously remarked, in this paper we were mostly interested in showing how double Higgs searches can be performed and interpreted in the general framework of EFT.We thus adopted a very simple analysis strategy avoiding unnecessary complications.It is clear, however, that the use of a more advanced procedure, as for instance a multivariate analysis, could be useful to improve the final sensitivity.In view of future high-energy colliders, it would also be interesting to investigate how to extract information about dimension-8 operators.A possible strategy, in this case, could be the use of angular distributions and a more systematic analysis of the high-energy tails of kinematic distributions.
Finally, a major point that is still missing in the literature is a full analysis combining the results obtained in the various double Higgs channels.Although γγb b is undoubtedly the easiest final state to look for, other channels with larger rates can be competitive and might lead to a higher sensitivity on the BSM coefficients.The estimates vary significantly in the literature depending on the details of the simulation and of the statistical analysis.The highest sensitivity is expected to come from the b bτ þ τ − final state, which has been claimed to lead to a determination of c 3 at the 40%-60% level at the highluminosity LHC with ffiffi ffi s p ¼ 14 TeV and L ¼ 3 ab −1 [5,9,13,54].A slightly worse sensitivity on c 3 has been estimated for the b bWW Ã channel [6,54], while extracting the trilinear coupling from the b bb b channel will be much harder.For the latter case it has been estimated that the signal can be distinguished from a background-only hypothesis with 95% probability if c 3 < 1.2 [11], as smaller values of c 3 lead to larger cross sections.For a fair comparison with our results, it should be noted that none of the above estimates, with the exception of those of Ref. [13], includes the uncertainty implied by the presence of other effective operators besides the one controlling the trilinear coupling.Our analysis suggests that this will have an important impact in all cases in which a high precision can be reached on the trilinear coupling.All analyses of double Higgs production, including the one illustrated in this work, agree in concluding that the observation of double Higgs production will be a challenge for the LHC, even at its high-luminosity upgrade.Extracting the Higgs trilinear coupling with sufficient precision will thus most certainly require combining as many different channels as possible.This will also open up the possibility of precisely testing NP by following the strategy outlined in this paper.Note added.-During the completion of this work Ref. [57] appeared, where an analysis of gg → hh → b bγγ is performed at a future 100 TeV collider in the context of the SM.The authors estimate a ∼40% sensitivity on the Higgs trilinear coupling with 3 ab −1 of integrated luminosity, which roughly agrees with our result.Their simulation differs in several aspects compared to ours: on the one hand, backgrounds with fake photons are included and their importance is pointed out; on the other hand, backgrounds are computed at LO (including the b bγγ continuum, whose simulation does not include extra jets at the ME level) and jet substructure techniques are not used.

APPENDIX A: ANGULAR DECOMPOSITION
By angular momentum conservation, the scattering amplitude can be decomposed into two contributions, M 0 and M 2 , mediating, respectively, J z ¼ 0 and J z ¼ AE2 transitions [58], The P μν 0;2 are (orthogonal) projectors where the form factors F □ , F Δ , and G □ are given in Ref. [59], and c gD0 , c gD2 are the coefficients of the higherderivative operators which give a next-to-leading correction to the nonlinear Lagrangian (4).They are related to the coefficients of the dimension-8 operators in (7) by the following simple relations: c gD0 ¼ ð4π=α 2 Þc gD0 , c gD2 ¼ ð4π=α 2 Þc gD2 .
APPENDIX B: SIMULATION

Cross sections of signal and backgrounds
The signal and background cross sections at the generation level are summarized in Table X.We simulate one extra jet from the matrix element in addition to the b bγγ system in all backgrounds except t th and γγjj.No generation-level cuts have been applied to the resonant backgrounds t th and Zh.The b bh samples are restricted to satisfy p T ðbÞ > 20 GeV at 14 TeV (increased to 30 GeV at 100 TeV) within jηðbÞj < 3 as well as ΔRðb; bÞ > 0.3.All extra partons are required to be within jηj < 5, and their p T thresholds are set to xqcut in MADGRAPH5.The size of γγjj is reduced by demanding p T ðjÞ > 15 GeV at 14 TeV (25 GeV at 100 TeV) within jηðjÞj < 5 and p T ðγÞ > 20 GeV (30 GeV at 100 TeV) within jηðγÞj < 2.5.Jets and photons are required to be separated in ΔR by 0.  approximation of the NLO one.This indicates that virtual corrections give a small contribution and the bulk of the NLO correction comes from real emissions.Another important aspect is the differential distributions.Although total rates agree well, the shape of distributions could differ between the NLO sample and the LO matched one, leading to different cut efficiencies.Figure 21 compares the distributions of several kinematic variables relevant for our analysis.While the higher value bins are subject to statistical fluctuations, overall the differential shapes of the two samples agree well.This justifies our use of the LO matched γγb b sample-whose generation requires much less CPU time-throughout our analysis.We expect that the same results will also apply in the case of a 100 TeV collider.

FIG. 3 (FIG. 4 (
FIG. 3 (color online).Differential cross sections obtained by including only the contribution of M 0 (dashed blue curve) or M 2 (dotted orange curve) in the SM, as functions of cos θ.Both curves are normalized to the total SM cross section.The partonic center-of-mass energy has been fixed to ffiffi ffi ŝ p ¼ 400 GeV in the left plot and to ffiffi ffi ŝ p ¼ 700 GeV in the right plot.

10 FIG. 5 .FIG. 6 .
FIG. 5. Left plot: Normalized differential cross section for pp → hh in the SM as a function of the invariant mass of the two Higgs bosons.The solid and dotted lines correspond, respectively, to ffiffi ffi s p ¼ 14 and 100 TeV.Right plot: Same as the left but with a logarithmic scale.

p
T> ðbÞ; p T> ðγÞ > 50 GeV; p T< ðbÞ; p T< ðγÞ > 30 GeV; 60 < m reco b b < 200 GeV; 60 < m reco γγ < 200 GeV; ð23Þ where p T> ðγÞ [p T< ðγÞ] denote the transverse momentum of the hardest [softest] photon, and p T> ðbÞ, p T< ðbÞ are similarly defined for the two b jets.At this stage, the broad mass windows in Eq. (23) are placed merely to allow for a fair comparison of the signal and background cross sections.

FIG. 7 (
FIG. 7 (color online).Angular distributions (in arbitrary units) after the cuts of Eq. (23) of the SM signal and the two dominant backgrounds, γγb b and t th, for ffiffi ffi s p ¼ 14 TeV.The dashed lines indicate the cuts of Eq. (24).

FIG. 8 .
FIG.8.Multiplicity of jets (left plot) and of reconstructed hadronic W bosons (right plot) after all cuts at ffiffi ffi s p ¼ 14 TeV.The solid, dotted, and dot-dashed curves denote, respectively, the SM signal, γγb b, and t th.

FIG. 9 .
FIG.9.Multiplicity of jets (left plot) and of reconstructed hadronic W bosons (right plot) after all cuts at ffiffi ffi s p ¼ 100 TeV.The solid, dotted, and dot-dashed curves denote, respectively, the SM signal, γγb b, and t th.

FIG. 10 .
FIG.10.Left: Pseudorapidity distribution of the most forward bottom quark or photon in the SM gg → hh → b bγγ process at 14 TeV (solid curve) and 100 TeV (dot-dashed curve).The vertical solid line corresponds to the acceptance region at the LHC.The dot-dashed vertical line denotes the acceptance region that should be adopted at 100 TeV to retain the same fraction of events of the 14 TeV case.Right: ΔR distance between two bottom quarks or photons in the SM gg → hh → b bγγ process at 100 TeV.The solid (dot-dashed) curve corresponds to the events with m hh > 1000 GeV (m hh < 1000 GeV).Both plots are made at the partonic level.

FIG. 11 .
FIG. 11.Number of signal events as a function of m recohh after all cuts at 14 TeV (left plot) and 100 TeV (right plot).The solid line is obtained by adopting the substructure technique with reduced photon isolation cone size (scenario Substructure II described in the text), whereas the dot-dashed line corresponds to a traditional jet-based analysis.

1 FIG. 15 (
FIG. 15 (color online).Breakdown of m hh categories in the plane ðc 2t ; c 3 Þ for the HL-LHC (left plot), and in the plane ðc 2t ; c 2g Þ for the FCC 100 (right plot).The various curves indicate the 68% probability contours for the following pairs of categories of Tables V and VII: 1 and 2 (dotted blue line); 3 and 4 (dot-dashed brown line); 5 and 6 (dashed red line).

22
suggestions and discussions.We would also like to acknowledge useful conversations with Jérémie Quevillon, Riccardo Rattazzi, Juan Rojo, and Andrea Wulzer.The work of R. C. was partly supported by the ERC Advanced Grant No. 267985 Electroweak Symmetry Breaking, Flavour and Dark Matter: One Solution for Three Mysteries (DaMeSyFla).The work of G. P. was partially supported by the Spanish Consolider-Ingenio 2010 Programme CPAN (CSD2007-00042), by CICYT-FEDER-FPA2011-25948, by the Severo Ochoa excellence program of MINECO under the grant SO-2012-0234 and by Secretaria d'Universitats i Recerca del Departament d'Economia i Coneixement de la Generalitat de Catalunya under Grant 2014 SGR 1450.

p 2 c p ν a p μ b p 2 T ðp a p b Þ − 2ðp b p c Þp ν a p μ c p 2 T ðp a p b Þ − 2ðp a p c Þp μ b p ν c p 2 T ðp a p b Þ þ 2p μ c p ν c p 2 T ; ðA2Þ and p 2 T
¼ 2ðp a p c Þðp b p c Þ=ðp a p b Þ − p 2 c is the transverse momentum of the Higgses.The expression of M 0 and M 2 is given by

10 FIG. 21 .
FIG. 21.Differential distributions for the b bγγ background obtained with the tree-level matched sample (solid lines) and the full NLO sample (dotted lines) for several variables at the 14 TeV LHC.Because of the mass window cuts at the generation level [see Eq. (B2)], the events are restricted to the mass window 60 < m reco b b < 200 GeV and 60 < m reco γγ < 200 GeV.
This estimate assumes that O gD0 and O gD2 are generated at 1-loop by the NP.The same assumption has been made on O g in the estimate of Eq. (3).
4=4, and O ¼ S 2 .One can thus identify m Ã ¼ m S and g Ã ¼ ffiffiffiffi ffi λ S p . 5 6R.C. would like to thank Riccardo Rattazzi for illuminating discussions on the validity of the effective theory in scattering processes.For a related analysis of WW scattering, see R. Rattazzi, contribution to "BSM physics opportunities at 100 TeV", CERN, 2014.

Table I
.

TABLE III .
Cut flow for the SM signal and the various backgrounds at ffiffi ffi s p ¼ 14 TeV.The values correspond to the number of events, assuming an integrated luminosity L ¼ 3 ab −1 .All reconstruction efficiencies and branching ratios to the final state under study are included.

TABLE IV .
Coefficients of the fit of the total signal rate [r ¼ σ × BRðhh → b bγγÞ] obtained after all cuts at 14 TeV and 100 TeV.The fit is based on a parametrization analog to that given in Eq. (16) for the cross section.The SM rates r SM include the NLO k factors of Eq. (17).

TABLE V .
Expected number of events after all cuts at ffiffi ffi s p ¼ 14 TeV in each of the six m reco hh categories considered, assuming an integrated luminosity L ¼ 3000 fb −1 .The last category is inclusive.
selection of events we increase the p T cuts on the b-tagged jets and photons while keeping the mass windows the same, p T> ðbÞ; p T> ðγÞ > 60 GeV; p T< ðbÞ; p T> ðγÞ > 40 GeV; 60 < m reco b b < 200 GeV; 60 < m reco γγ < 200 GeV:

TABLE VII .
Expected number of events after all cuts at ffiffi ffi s p ¼ 100 TeV in each of the six m reco hh categories considered, assuming an integrated luminosity L ¼ 3000 fb −1 .The last category is inclusive.The numbers in the parentheses are obtained by removing the veto on hadronic W's.

TABLE VIII .
Comparison between the traditional jet-based analysis (second line) and the substructure analysis (third, fourth, and fifth lines) on SM signal events.The three scenarios, Substructures I, II, and III, differ by the treatment of the γγ system and are defined in the text.

TABLE IX .
Expected numbers of signal and background events at ffiffi ffi s p ¼ 100 TeV after the jet substructure analysis (scenario Substructure II described in the text with photon isolation cone R iso ¼ 0.2) assuming an integrated luminosity L ¼ 3000 fb −1 .The last category is inclusive.
FIG.19(color online).Probability distributions as functions of c6 at the HL-LHC.The black continuous curve refers to the standard exclusive analysis discussed in Sec.IV.The dotted red curve refers to an inclusive analysis in the left plot, and to an analysis without marginalization on cH , cd , cg , and cu in the right plot.See text.FIG.20 (color online).Width of the 68% probability interval on c6 (orange region) as a function of the uncertainty on cu for the FCC 100 .See text.

TABLE X .
Cross sections (in units of fb) for signal and backgrounds at 14 TeV and 100 TeV after generation-level cuts (see text for details).