I Introduction

CERN-PH-TH-2015-015

Effective field theory analysis of double Higgs

[0.15cm] production via gluon fusion

Aleksandr Azatov, Roberto Contino, Giuliano Panico and Minho Son

Theory Division, Physics Department, CERN, Geneva, Switzerland

Institut de Théorie des Phénomènes Physiques, EPFL, Lausanne, Switzerland

IFAE, Universitat Autònoma de Barcelona, Bellaterra, Barcelona, Spain

Abstract

We perform a detailed study of double Higgs production via gluon fusion in the Effective Field Theory (EFT) framework where effects from new physics 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 channel and perform a thorough simulation of signal and background at the TeV LHC and a future 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 accuracy can be reached on the trilinear coupling at a future TeV collider with . Only an determination seems instead 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 GeV and properties similar to those predicted for the Standard Model (SM) Higgs boson has been a major leap forward in this direction Aad et al. (2012); Chatrchyan et al. (2012). 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 New Physics (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 . They set limits on possible modifications of the Higgs couplings, whose naive expected size is of order , where is the mass of the new states, is their coupling strength to the Higgs boson and a SM coupling. The improved performances in energy and luminosity which will characterize the LHC Run2, on the other hand, give the opportunity to directly probe the EWSB dynamics at much higher energies () through the study of scattering processes. For a typical scattering energy , effects from New Physics are expected to be of order , hence enhanced by a factor 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 inferring how much theoretical space is being probed.

Double Higgs production via gluon fusion is one example of scattering process which can disclose key information on the EWSB dynamics, its underlying symmetries and 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 Baur et al. (2003, 2004); Dolan et al. (2012); Papaefstathiou et al. (2013); Baglio et al. (2013); Yao (2013); Barr et al. (2014a); Barger et al. (2014a); Ferreira de Lima et al. (2014). 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 TeV LHC and a future TeV proton-proton collider by focusing on the final state. We perform a detailed Montecarlo (MC) simulation of signal and background, and use the kinematic distributions to maximize the sensitivity on different effective operators. Previous studies in the EFT context appeared in Refs. Contino et al. (2012); Goertz et al. (2014), and a first analysis of the impact of New Physics on the kinematic distributions can be found in Ref. Chen and Low (2014).

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 TeV collider, giving a first estimate of the sensitivity on the kinematic tail at large invariant masses of the Higgs pair. Our Montecarlo analysis of the decay mode is illustrated in Section IV, while results on the sensitivity to the coefficients of the effective Lagrangian are collected in Section V. We compare with previous analyses of the final state in Section 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 background into the Appendices.

## 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 which 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 is linearly realized at high energies. It thus referred to as the “linear” Lagrangian. In the second, more general formulation, is non-linearly realized, hence the name of “non-linear” 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. Contino et al. (2013), 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 non-linear 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:

 Llin=LSM+ΔL6+ΔL8+… (1)

The lowest-order terms coincide with the usual SM Lagrangian , whereas contains the deformations due to operators of dimension , with 1 For our purposes it is sufficient to focus on the operators involving the Higgs boson. The relevant ones in are 2

 ΔL6⊃¯¯cH2v2∂μ(H†H)∂μ(H†H)+¯¯cuv2yt(H†H¯¯¯qLHctR+h.c.)−¯¯c6v2m2h2v2(H†H)3+¯¯cgg2sm2WH†HGaμνGaμν, (2)

where GeV and GeV is the physical Higgs mass. In order to classify the various operators in our effective Lagrangian and estimate their coefficients we adopt the SILH power counting Giudice et al. (2007). This is based on the assumption that the NP dynamics can be broadly characterized by one mass scale , at which new states appear, and by one coupling strength . 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 , while each additional covariant derivative is suppressed by . If the Higgs is a composite Nambu-Goldstone (NG) boson Kaplan and Georgi (1984), the operators , and can be generated only through some small explicit breaking of the global invariance since they violate the Higgs shift symmetry Giudice et al. (2007)3 In this case the SILH estimate of the coefficients in Eq.(2) is:

 ¯cH,¯cu,¯c6∼v2f2,¯cg(4πα2)∼v2f2×λ2g2∗, (3)

where is a weak spurion coupling suppressing , while the suppression of and is controlled respectively by and . In realistic models, the minimum amount of explicit breaking entering is given by the top Yukawa coupling , so that one expects . For example, one has in models with fully composite , whereas partial compositeness of both and leads to  Giudice et al. (2007).

In the case of the non-linear Lagrangian, the terms relevant for our analysis are

 Lnon−lin⊃−mt¯tt(cthv+c2th2v2)−c3m2h2vh3+g2s4π2(cghv+c2gh22v2)GaμνGaμν, (4)

where denotes the physical Higgs field (defined to have vanishing vacuum expectation value). Compared to the linear Lagrangian, the couplings of Eq. (4) effectively resum all the corrections of order . Therefore, the non-linear Lagrangian only relies on the derivative expansion, in which higher-order terms are suppressed by additional powers of . The linear formulation further requires 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:

 ct=1−¯¯cH2−¯¯cu,c2t=−12(¯cH+3¯cu),c3=1−32¯¯cH+¯¯c6,cg=c2g=¯¯cg(4πα2), (5)

where . Notice in particular that the same operator which gives a (non-universal) modification of the top Yukawa coupling also generates the new quartic interaction .

It is worth at this point to comment about the normalization of the operator in Eq. (2). This differs from the one of Ref. Contino et al. (2013), where it was defined (here we introduce the symbol to distinguish it from in Eq. (2)), with denoting the coefficient of the operator in the SM Lagrangian. We find , which should be compared with the relation between and (valid at all orders in ) appearing in Eq. (5). For small values of one has , hence , so that the dependence of upon is the same as in Eq. (5) up to small corrections. The parametrization of Eq. (2) is thus more convenient than that of Ref. Contino et al. (2013) when (or ) becomes large, since the formula for the trilinear coupling remains linear in . This is relevant for this analysis, since we anticipate that the results of Section V will constrain values of larger than 1 at the LHC.

### ii.1 Modified power counting for the Higgs trilinear coupling

The estimates of Eq. (3) shows that assuming the SILH power counting the modification of the Higgs trilinear coupling is expected to be of order , 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 boson– of some new strong dynamics. In this case no weak spurion suppression is required to generate , and the naive estimate of is enhanced by a factor compared to the SILH case, where we conveniently defined . The Higgs trilinear coupling (in SM units) is thus expected to be of order , and large modifications are possible for large even if . For example, gives and , corresponding to . 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 . This is in addition to the tuning which must occur in the vacuum alignment even for a NG boson Higgs. Notice that similarly to , other operators which previously required a breaking of the Goldstone symmetry to be generated will now be unsuppressed. In particular, and the coefficient of the operator get enhanced by a factor compared to their SILH estimate, thus leading to modifications of the Higgs couplings to gluons and photons of order times their (loop-induced) SM value. However, for the current LHC constraints on these couplings can be easily 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  and a coupling strength , couples to the SM sector only through the Higgs mass term portal: . The Higgs field is thus elementary, and couples to the composite operator , made of fields of the strong sector, with coupling . This leads to the following estimates:

 Δλ4∼λ2g2∗,¯cH∼λ2g4∗v2f2,¯cu∼λ2g4∗v2f2y2t16π2,¯c6∼λ3g4∗¯λ4v2f2, (6)

where is the correction to the coefficient of the operator induced by the strong dynamics. The contribution to follows from a Higgs-dependent modification of the top quark kinetic term arising at the 1-loop level, and it is subleading compared to . 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 . An additional factor should be included in the above estimates if they arise from diagrams involving a loop of strongly-coupled particles. 4 Since , one can have an enhanced NP correction to the Higgs trilinear coupling for . However, requiring that no fine-tuning occurs in the coefficient of the quartic operator –i.e. in the physical Higgs mass– implies . Furthermore, only for 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 . In this scenario it is thus possible to avoid tuning the Higgs mass but the corrections to the Higgs trilinear are at most of ; larger enhancements require fine-tuning. For example, , GeV and gives and .

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, can be of order 1 or larger without invalidating the perturbative expansion in powers of the Higgs field.

### ii.2 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 . This means that if the linear Lagrangian cannot be consistently used anymore, although we can still rely on the non-linear one. As regards the derivative expansion, the same considerations apply instead to both descriptions. The perturbative parameter in that case is , so that naively higher-order operators become important at , 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 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:

 (7)

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: 5

 ¯cgD0,2(4πα2)∼v2f2×m2Wm2∗. (8)

We can now compare the contributions of the various operators to double Higgs production via gluon fusion. 6 By using the estimates of Eqs. (3),(8), the scattering amplitude can be schematically expressed as

 A(gg→hh)∼(αs4π)×[y2t(1+O(v2f2))+g26(E)+g28(E)+…], (9)

where

 g26(E)∼¯cg4πα2E2v2∼λ2E2m2∗,g28(E)∼¯cgD0,24πα2E4v2m2W∼g2∗E4m4∗. (10)

The first term in the square parentheses of Eq. (9) corresponds to the SM contribution plus the correction implied by and 7 The second and third terms in the parentheses denote the contributions of respectively and , , and thus define the coupling strengths of the interactions mediated by these operators. Since these couplings grow with , 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 , hence always smaller than the SM term. We thus see the potential advantage of studying 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 New Physics dynamics in the regime in which it gives large effects Contino et al. (2014). In practice, as we will discuss more in detail in the following, such regime is difficult to exploit in the process considered in this work. This is due to the steep fall off of the gluon pdfs at large and to the small final-state branching ratio, which strongly limit the exploration of the region with large invariant mass. For the coupling is the first to become larger than at a scale . The smaller is, the higher the crossover scale becomes, until is reached where . Above this energy, for , the largest effects come from the dimension-8 operators. The validity of the effective Lagrangian extends up to energies , so that the coupling strengths are bounded by

 g6(E)<λ≲g∗,g8(E)

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 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 for ; at this scale . Hence, a further condition to be satisfied in an analysis which includes only dimension-6 operators neglecting those with dimension 8 is:

 g6(E)<λ2g∗. (12)

Let us then indicate with the precision obtained in such analysis on , by making use of events with invariant masses up to . This means that the smallest value of probed is , so that Eqs. (11),(12) can be re-cast into the following constraints:

 λ >gmin (13) λ2g∗ >gmin. (14)

These inequalities define the region in the 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 , obtained by increasing the precision for a fixed energy , lead to a larger viable region. Notice that the two cases (fully composite ) and (partially composite and ) can be probed only if . This is not the case of course if the analysis does not have enough sensitivity to probe the SM cross section.

### ii.3 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 non-linear 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 be then 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 region.

The Feynman diagrams that contribute to the process are shown in Fig. 2.

Each diagram is characterized by a different scaling at large energies . We find:

 A□∼c2tαs4πy2t,A△∼ctc3αs4πy2tm2h^s(logm2t^s+iπ)2,A△nl∼c2tαs4πy2t(logm2t^s+iπ)2,A3∼cgc3αs4πm2hv2,A4∼c2gαs4π^sv2, (15)

where , are the amplitudes of respectively the box and triangle diagram with Higgs exchange, is the amplitude of the triangle diagram with the interaction, and and are the amplitudes of the diagrams with the Higgs-gluon contact interactions. One can see that each NP contribution affects the distribution in a different way. In particular, the diagrams that depend on the Higgs trilinear coupling are always suppressed at large , and their contribution affects the process mostly at threshold. Modified values of the top Yukawa coupling and the non-linear interactions and , instead, tend to increase the cross section at higher invariant masses. Finally, including the dimension-8 operators would lead to an additional contribution to growing as and distort the tail of the 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 as the main kinematic variable to characterize signal events.

By focussing on gluon fusion, the total cross section for the process can be written as a simple polynomial of the parameters of the effective Lagrangian: 8

 σ=σSM[A1c4t+A2c22t+A3c2tc23+A4c2gc23+A5c22g+A6c2tc2t+A7c3tc3+A8c2tctc3+A9c2tcgc3+A10c2tc2g+A11c2tcgc3+A12c2tc2g+A13ctc23cg+A14ctc3c2g+A15cgc3c2g]. (16)

The LO value of the numerical coefficients and of the SM cross section is reported in Table 1 for hadron colliders with center-of-mass energy TeV and TeV. 9 They were computed with our dedicated C++ code linked to QCDLoop Ellis and Zanderighi (2008) and to the LHAPDF routines Whalley et al. (2005), by setting GeV, GeV and using the CTEQ6ll parton distribution functions. The factorization and renormalization scales have been fixed to .

Besides , the other kinematic variable that characterizes the 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, and , describing transitions with respectively and , where is the projection of along the beam axis. The amplitude receives a contribution only from the box diagram and from the dimension-8 operator (through the last diagram of Fig. 2). All the other diagrams, instead, mediate transitions. The explicit expression of and is reported in Appendix A for convenience. While transitions with all integer values of the total angular momentum can occur, the leading contributions to and come respectively from (s-wave) and (d-wave). A simple angular momentum decomposition then shows that, up to small corrections, the angular dependence is , 10 In the SM, the contribution of to the cross section is always small, as illustrated by Fig. 3: it is negligible at the peak of the distribution (GeV) and smaller than on its tail (GeV).

A shift in the top Yukawa coupling due to New Physics 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 can give a sizable contribution to through the last diagram of Fig. 2. If this were the case, one could reach a better sensitivity on 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 in the plane . We defined to be the ratio of the cross section induced by alone, , to the total cross section, , obtained by adding the contribution of to the SM. Both cross sections are computed at the partonic level for the process and by integrating over angles . We set the coefficient equal to its NDA estimate (8) and choose as benchmark values GeV, TeV, which correspond to , .

For the contribution of is as large as the SM one. The plot of Fig. 4 shows that this occurs at the crossover scale TeV, in agreement with the naive estimate of section II.2, . Increasing the value of (i.e. decreasing ) enhances the importance of the component of the cross section, hence that of , but the effect is never large. For example, the value of the crossover scale is reduced at most by . We thus conclude that, although an analysis of angular distributions could in principle help disentangling the effects of dimension-8 operators, in practice there is little leverage, and the efficacy of such strategy is further reduced in the case of the final state by the limited range in which can be realistically probed. In our analysis of we will thus make use of the distribution as the main tool to probe the effects induced by New Physics, neglecting for simplicity angular distributions.

## Iii Double Higgs phenomenology at 14TeV and 100TeV

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 and . 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 TeV LHC, where the NLO total SM production rate is around , but remains true even at a hypothetical future TeV machine, where the enhanced rates are unavoidably accompanied by larger backgrounds. In order to obtain a large enough branching ratio, at least one of the two Higgses must decay into a pair. For the second Higgs different choices seem possible and have been considered in the literature, namely , , and . The channel has the highest signal rate ( in the SM), but is plagued by a large QCD background. Even imposing four -tags, it seems hard to exploit and could require a sophisticated search strategy Baur et al. (2003); Dolan et al. (2012); Ferreira de Lima et al. (2014). The channel with the second highest rate is , with a branching ratio in the SM. Its observation is also threatened by a large background, mainly coming from  Dolan et al. (2012); Baglio et al. (2013). It has been recently claimed that a good signal-to-background 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 Papaefstathiou et al. (2013). Although encouraging, the results of this analysis suggest that the final state can be observed at the TeV LHC only with its high-luminosity extension. The other channel which has been extensively studied in the literature is  Baur et al. (2003); Dolan et al. (2012); Barr et al. (2014a); Baglio et al. (2013); Goertz et al. (2014). It is very promising and potentially relevant for the TeV LHC, since its SM branching ratio is sizable, , 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 will have to be fully assessed by an experimental study.

For the purpose of our study we will focus on the channel ( in the SM), which has been considered to be the cleanest one despite its small rate. As shown by previous theoretical studies Baur et al. (2004); Barger et al. (2014a); Yao (2013); Baglio et al. (2013) as well as a recent non-resonant experimental search at TeV Aad et al. (2014a) and a study for the TeV high-luminosity LHC  ATLAS Collaboration (2014) 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 effective field theory 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 TeV LHC and a future TeV collider.

It is well known (see for example Ref. Contino et al. (2012)) 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 (). The peak of the distribution is essentially determined by the fast decrease of the gluon parton distribution functions and is located around . This conclusion is valid independently of the collider center-of-mass energy . 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 TeV machine is around times bigger than the one at the TeV LHC), with small effects on the distribution. As it can be seen from Fig. 5, this latter is modified significantly only on its tail (), which is enhanced at larger collider energies due to the higher luminosity of gluons. In fact, as we will discuss in section IV.1, an important change is also present in the pseudo-rapidity distribution of the Higgs bosons, which is related to the amount of boost determined by the initial parton energies.

Due to the above considerations, one would naively expect that an analysis similar to the one designed for TeV would continue to work well at 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, due to the larger cross section and enhanced tail, much higher values in 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 TeV collider requires some modification of our TeV analysis strategy in order to better reconstruct events with a higher boost. Although we will not present a fully optimized analysis for the 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 and 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 new physics contributions to the cross section. To get a rough approximation, however, we can estimate the reach by demanding that a few events (we choose for the estimate) are still present in the tails of the SM distributions above the considered (or ) value. To be more realistic we also assume a overall signal efficiency due to kinematic cuts. The estimates can be easily extracted from the plots of Fig. 6, which show the integrated differential cross sections for the SM signal as a function of the lower cut on and of the minimum .

By including the branching ratio of the most important decay modes we obtain the results shown in Table 2.

The reach on is rather limited in the channel due to the small signal rate, and even at a TeV collider it does not extend much beyond 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 . The maximal value of 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 . In order to resolve the jets with the standard techniques we must demand , where and are the reconstruction cones used by the ATLAS and CMS collaborations respectively. It is easy to see that for 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 2 show that the size of the boosted phase space at TeV is almost negligible for the channel and very limited for , and . The story changes however for an TeV collider. In this case the SM signal can be probed up to in the channel and up to for . 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 details in Section IV.1.

## Iv Analysis of pp→hh→γγb¯b

After discussing the general properties of the process, we now focus on the final state and describe our analysis strategy. We generated the signal events through our own dedicated code described in section II.3. 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 PDF’s (LO PDF with LO ), setting the factorization and renormalization scales to and GeV, GeV, we find the following LO cross sections for the Standard Model: fb and fb for respectively TeV and TeV collider energies. In order to (partially) include the NLO and NNLO corrections, we rescale the cross sections by the k-factors

 k14TeV=2.27,k100TeV=1.75 (17)

which were computed for the SM in Ref. de Florian and Mazzitelli (2013)12

The main backgrounds that we considered are the non-resonant processes and and the resonant processes , and (with the subsequent decays , , and ). 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_aMCNLO v2.1.1 Alwall et al. (2014) by switching off virtual corrections (i.e. working in LO mode). The output has been interfaced with PYTHIA v6.4 Sjostrand et al. (2006) for parton showering and hadronization, and with FastJet v3.0.6 Cacciari et al. (2012) 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 , and have been matched up to one additional jet via the -jet MLM matching Alwall et al. (2008) to partially account for NLO effects. In the case of we found that allowing for the extra jet increases the cross section by a factor . We found that this is in good agreement with a complete NLO computation performed with MadGraph5_aMCNLO, indicating that real emissions in this case represent the bulk of the NLO correction. Considering that 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 and a thorough discussion of the NLO correction are given in Appendix B.2. For the resonant background we apply instead a simple k-factor rescaling of the LO cross section to the NLO value reported in Ref. Dittmaier et al. (2011). Finally, the background was simulated without matching and a k-factor was applied (this is a more conservative choice than the one proposed in Alwall et al. (2014) 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 ratio is approximately the same for both the SM and in a generic point of the BSM parameter space. The rescaling is performed individually for each of the bins in of Tables 57 and 8. Our analysis takes correctly into account possible non-universal 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 Section IV.1, 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 multiply 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 , 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 level in the high- categories used for the boosted analysis.

We designed different strategies for TeV and TeV colliders. We begin by describing the one at TeV. Events are triggered by demanding exactly two isolated photons satisfying the minimal reconstruction requirements

 pT(γ)>25 GeV,|η(γ)|<2.5. (20)

To forbid extraneous leptonic activity we veto events with isolated leptons (electrons or muons) satisfying

 pT(l)>20 GeV,|η(l)|<2.5. (21)

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 (0.3) satisfies (0.85). We have checked that more sophisticated photon isolation criteria, like the one proposed in Frixione (1998), give similar results. We then cluster the events into anti- jets Cacciari et al. (2008). We accept only jets with

 pT(j)>25 GeV,|η(j)|<2.5 (22)

and require that at least two of them are -tagged (the leading two -jets are selected if more than two -jets exist). We assume 70% efficiency for tagging , corresponding to 1% of mis-tag rate  Chatrchyan et al. (2013); CMS Collaboration (2013a); ATLAS Collaboration (2013a)13 and 80% efficiency for photon tagging  CMS Collaboration (2013b); ATLAS Collaboration (2013a) 14. After this step, an event consists of two isolated photons, two -tagged jets, and possible additional jets (whether or not -tagged).

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 -tagged jets satisfying

 pT>(b),pT>(γ)>50 GeV,pT<(b),pT<(γ)>30 GeV,60

where () denote the transverse momentum of the hardest (softest) photon, and , are similarly defined for the two -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.

In the signal, the and subsystems tend to be in a back-to-back configuration, and the angular distance between two photons (and similarly between two -jets) is of order 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 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 , while retaining most of the signal. 15 As a final cut, we restrict the invariant masses of two -jets and of the two photons to the Higgs mass window 16

 105

The resulting cut flow is shown in Table 3.

We report in Table 4 a fit of the signal rate () obtained after all cuts based on a parametrization analog to that given in Eq. (16) for the cross section.

With our analysis strategy, and are the two major backgrounds. The latter tends to produce extra jets from the hadronic decay of the ’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 , in addition to all the previous cuts, can remove roughly 80% of the background while keeping 70% of the signal. Alternatively, one could look for hadronic ’s by iteratively forming jet pairs with invariant mass lying in a given window around the mass, for example we will use GeV in the following. As illustrated in Fig. 8, vetoing hadronic ’s in addition to our cuts can reduce 50% of the 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 TeV, motivated by the necessity of retaining as many signal events as possible.

As discussed in Section II.3, several diagrams with different energy scalings contribute to the signal. We can thus use the invariant mass of the reconstructed 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 . The corresponding numbers are reported in Table 5 and will be used in section V to extract our bounds on the coefficients of the effective operators.