Higgs boson pair production in the D=6 extension of the SM

We derive the constraints that can be imposed on the dimension-6 effective theory extension of the Standard Model, using gluon fusion-initiated Higgs boson pair production at the LHC. We use a realistic analysis focussing on the $hh \rightarrow (b\bar{b}) ( \tau^+ \tau^- )$ final state, including initial-state radiation and non-perturbative effects. We include the statistical uncertainties on the signal rates as well as conservative estimates of the theoretical uncertainties. We first consider a theory containing only modifications of the trilinear coupling, through a $c_6 \lambda\, H^6/ v^2$ Lagrangian term, and then examine the full parameter space of the effective theory, incorporating current bounds obtained through single Higgs boson measurements. We also consider an alternative scenario, where we vary a smaller sub-set of parameters. Allowing, finally, the values of the other coefficients to vary within \textit{projected} experimental ranges, we find that the currently unbounded parameter, $c_6$, could be constrained to lie within $|c_6| \lesssim 0.6$ at 1$\sigma$ confidence, at the end of the high-luminosity run of the LHC (14~TeV) in the full model, and to $-0.6 \lesssim c_6 \lesssim 0.5$ in the alternative model. This study constitutes a first step towards the inclusion of multi-Higgs boson production in a full fit to the dimension-6 effective field theory framework.


Introduction
The scalar particle recently discovered by the ATLAS [1] and CMS [2] collaborations at the Large Hadron Collider (LHC) appears to be compatible with the Higgs boson of the Standard Model (SM) [3][4][5]. In particular, it seems to behave like a CP-even scalar, with couplings to the gauge bosons and fermions that agree at the O(20-100)% level [6-8] with those predicted by the SM, in the cases where the experiments are already sufficiently sensitive (gauge bosons and third generation fermions). The couplings to the SM fields have been probed via the production and decay modes of the Higgs boson and this is a programme that will continue in future runs of the LHC and forthcoming colliders.
The story is different, however, in the "pure Higgs" sector, characterised by the following (D ≤ 4) potential post-electroweak symmetry breaking (EWSB): where the self-couplings λ =λ = m 2 h /2v 2 within the SM, with v 246 GeV the vacuum expectation value. 1 Specifically, only the first term of Eq. (1.1) has been probed, through the measurement of the Higgs boson mass, m h 125 GeV. Hence, direct determination of the terms proportional to h 3 and h 4 is an essential experimental measurement that will provide access to new phenomena, such as a richer scalar sector or heavier coloured particles, or (in a rather grim scenario) put the validity of the SM on even more solid grounds.
At colliders, terms proportional to h n can be probed through the simultaneous production of (n−1) Higgs bosons. 2 Unfortunately, the production rates for these processes are not as big as for the single Higgs boson production processes. This is due to the relatively large invariant mass of the final state system and the fact that the processes are loop-induced. At the LHC with 14 TeV proton-proton centre-of-mass energy, triple production is expected to be rather rare, with a cross section of O(0.1 fb). This renders any measurement of the coefficient of the h 4 term hard, if not impossible, even at the high-luminosity LHC (HL-LHC) [10][11][12]. The prospects for Higgs boson pair production are significantly better, with a SM cross section almost three orders of magnitude larger (about 30-40 fb [13][14][15][16][17][18][19]). However, this process is still particularly challenging to detect, both in the SM [20][21][22][23][24][25][26][27][28][29][30][31] and beyond [9,. Moreover, sensitivities to the double Higgs boson production process do not immediately translate to sensitivities to the trilinear coupling. The reason is that several diagrams contribute to the production process, but only few of them involve the coefficients λ, that are furthermore associated with an off-shell Higgs propagator.
Since no new particles beyond the SM have been observed by the LHC experiments so far, they are either well hidden, weakly coupled, or simply heavier than a (couple of) TeV. In the latter case one can adopt an effective field theory (EFT) approach, where the effects of the high-scale physics are parameterized by a set of higher-dimensional operators, suppressed by a large mass scale Λ. In the present paper we investigate the Higgs boson pair production (hh) in the framework of the dimension-6 effective field theory (D=6 EFT) extension of the SM, modifying in particular (but not only) the Higgs potential (1.1). This allows us to derive model-independent constraints on new physics that may be beyond direct reach. The hh process can be a source of additional meaningful information, cutting through regions of the parameter space of the EFT coefficients in non-trivial ways. For example, it may help clarify whether or not the Higgs boson is really part of an SU (2) L doublet. 3 Studying the hh process is essential to probe the full set of the coefficients in the D=6 EFT extension of the SM. In particular, the coefficient of the D=6 operator H 6 currently remains unconstrained by single Higgs boson measurements, but can contribute to the hh process through direct modification of the Higgs boson trilinear self-coupling. This paper is organised as follows: in section 2 we examine the EFT Lagrangian containing a complete set of relevant dimension-6 operators that we wish to investigate. In section 3 we focus on the terms relevant to gluon-fusion-initiated Higgs boson pair production after EWSB and relate them to the SM EFT, i.e. the SM with the top quark integrated out. In section 4 we examine the impact of the dimension-6 operators on the decays of the Higgs boson and in section 5 we present our set-up for the analysis of the key process pp → hh → (bb)(τ + τ − ), that we then employ explicitly as an example of our framework to generate constraints. We conclude in section 6. We provide additional information on our conventions in appendix A. Appendix B provides technical details on the derivation of the Lagrangian after EWSB.

Higgs boson effective theory
New Physics associated to a new scale Λ v can be described in a model-independent way by augmenting the Lagrangian of the SM with all possible gauge-invariant operators of mass dimension D > 4, where the leading effects arise from D = 6 operators (neglecting lepton-number violating operators, irrelevant to our study). Working at this level, the extension of the SM that we will consider for our analysis of Higgs boson pair production reads where α s is the strong coupling constant and α ≡ g 2 /4π. Compared to the full set of D = 6 operators that can be formed out of the SM field content, as first obtained in [56] and reduced to a non-redundant minimal set in [57], we neglect a class of operators whose effect is already constrained from precision tests to be at most 1% with respect to the SM [58,59] (see also [60]). Including these operators would have a negligible numerical impact on the analysis, given the experimental and theoretical errors. 4 Precision measurements also lead to the approximate restrictions [58] which we will employ in the following. Thus our setup corresponds to a restricted SILH Lagrangian [61] where c T has been set to zero and the relations (2.2) are used. We also assume minimal flavour violation (MFV) [62], which leads to the coefficients of the Yukawa-like terms in the second row of Eq. (2.1) being proportional to the (SM-like) Yukawa couplings y t,b,τ and in particular allows to neglect the corresponding contributions involving the light fermions, hence Q L = (t L , b L ) and L L = (ν τ , τ L ). Note that the latter can also be justified without assuming MFV [63]. In particular, this helps to exclude the possibility of largely modified hierarchies in fermion-Higgs couplings with respect to the SM, which might have also spoiled the hierarchies between the different Higgs-production mechanisms. In order to further simplify the parameter space, we will also set c τ = c b for the rest of the paper. 5 In the spirit of the SILH model we have normalised the operators c g,γ including a loop suppression factor: in any perturbatively-decoupling extension of the SM these operators can only be generated at the loop level. Moreover, this is also convenient when comparing the D = 6 theory with the top quark EFT, that is, the limit where the top quark is integrated out from the SM, as we will discuss in section 3. After EWSB several operators might contribute to the same interaction and field redefinitions need to be introduced in order to have canonically normalised kinetic terms. We examine the terms relevant to Higgs physics in the next section.

Relevant Lagrangian terms
The dimension-6, CP-even, operators listed in the Lagrangian of Eq. (2.1) that affect the production of multiple Higgs bosons via gluon fusion (at leading order) are: where in the first line we have included the corresponding SM operators that will receive additional contributions in the effective theory. 6 We now write where T represents the 3 generators of SU (2) L , ξ represents the 3 Goldstone degrees of freedom that will be absorbed by the gauge bosons, and h is the physical Higgs boson. To have canonical normalisation of the Higgs field, we choose to perform the field redefinition 7 We further redefine c i → c i Λ 2 /v 2 to absorb the suppression factor into the c i coefficients.
We thus obtain the following interactions in terms of the Higgs boson scalar h, relevant to Higgs boson pair production: where ∆ = c 6 + 3c H /2. From the above, it can be seen that the SM relation of λ =λ is broken by the EFT effects: an accurate measurement of both couplings is thus a powerful probe of new physics in the Higgs sector, although, as already mentioned, measurement of the quartic coupling does not seem to be possible in the foreseeable future.

From SM EFT to dimension-6 EFT
It is useful to compare and contrast the dimension-6 extension of the SM with the EFT that results from taking the top mass to infinity within the SM framework. This will help us in writing down the cross section formula for gg → hh in the D=6 EFT. There are several modifications necessary to incorporate the effect of the D=6 EFT operators in Higgs boson pair production via gluon fusion (see Fig. 1): • The top and bottom quark Yukawa couplings will be modified according to the third line in Eq. (3.4). These modifications appear in diagrams 1A and 1C.
• The new four-point fermionic interactionsf f h 2 will introduce a new 'triangle' diagram, with two Higgs bosons produced at the apex of the triangle. These are described by the fourth line in Eq. (3.4) and appear in diagram 1E.
• Two new effective theory "tree-level" diagrams will contribute since we now have additional gluon-gluon-Higgs (diagram 1B) and gluon-gluon-Higgs-Higgs interactions (diagram 1D), according to the second line in Eq. (3.4).
These are relatively straightforward to incorporate into a Monte Carlo event generator. The shift in the Higgs boson couplings to fermions and to gauge bosons will also lead to changes in its branching ratios. We discuss those changes in section 4. Our goal is to modify the Standard Model matrix element with finite top mass to readily incorporate the modifications coming from the EFT. We first note that in the heavy top mass limit, the 'low energy theorem' states that the interactions of one or two Higgs bosons with two gluons are given by: It is interesting to point out that the gghh interaction that remains in the heavy top mass approximation corresponds to a spin-0 di-gluon state in gg → hh, while the spin-2 contributions vanish in this limit. In the effective theory language, the spin-2 contributions correspond to operators of higher dimensionality. The differential partonic cross section for gluon-fusion-initiated Higgs boson pair production in the SM, with the full top mass dependence, is given by [16]: where, in the SM, F , F 2 and G 2 are form factors, given e.g. in Ref. [16], in which G 2 corresponds to a spin-2 contribution, whileŝ andt are the usual Mandelstam invariants. In the limit of a large quark mass in the loop, m Q ŝ, the form factors reduce to: Note that the function G 2 is sub-dominant in this limit, in correspondence with the fact that the spin-2 terms are absent in Eq. (3.6). We now derive, starting from Eq. (3.7), the cross section for the hh process in the D=6 EFT. The complete set of diagrams is shown in Fig. 1. Using the above limiting values of the form factors, one can re-write the Lagrangian of Eq. (3.6) as: where F hq = −F hq 2 = 2/3 are the values of the form factors in the heavy quark effective theory: We begin by considering the correspondence between diagrams 1A and 1B as well as between 1C and 1D. This corresponds to comparing equivalent terms in Eq. (3.10) and Eq. (3.4). We can immediately conclude that the following identifications can be made at the amplitude level to obtain the contributions of the pure EFT diagrams: To obtain diagram 1E, one needs to essentially 'remove' the propagator from a diagram of type 1A, while keeping the dependence on the quark mass in the triangle loop via the full form factor F . This can be done by multiplying the factor which includes a combinatoric factor of 2, f = t, b, and (3.14) The necessary additional modifications correspond to trivial replacements of the triple-Higgs coupling and the Yukawa coupling in the SM-like diagrams, according to Eq. (3.4). After all these modifications, we now arrive at the parton-level differential cross section for the process gg → hh in the D=6 EFT: where we have suppressed the bottom quark contributions for simplicity.

Higgs boson decays in dimension-6 EFT
We now move forward to study the impact of the operators in Eq. (2.1) on the decays of the Higgs bosons. In Table 1 we provide an overview on which coefficients enter the various decays at the tree level (second column), at the one-loop level, considering only QCD corrections to the insertions (third column), as well as at the full one-loop level (fourth column). Here, we focus on the key decays h → bb, h → τ τ , h → W W , and h → γγ, that are in particular important for the analysis of Higgs boson pair production. In our numerical study we will, however, consider all significant decays in the EFT. We note that even though operators may not enter a given pp → hh → (xx)(yy) final state, they will still be relevant since they change the overall branching ratios. We also include in Table 1 the coefficients entering gg → h and gg → hh production for completeness. It is interesting to observe that several operators may enter both production and decay, introducing non-trivial correlations for the behaviour of the cross section of a given final state mediated through Higgs boson pair production. Finally, note that the presence of the coupling c γ can lift the h → γγ decay from one-loop in the SM formally to tree-level in the EFT. The same is true for c g concerning gg → h and gg → hh. On the other hand, in any perturbatively-decoupling extension of the SM the operators c γ , c g , c HW and c HB can only be generated at the loop level. Thus, we will in particular not insert them into loop diagrams at the order considered [68].
In the present article we employ the eHDECAY code [69] to calculate the branching ratios of the Higgs boson according to our EFT formalism. The program eHDECAY contains the leading-order SM contributions, the NLO QCD contributions, the NLO EW corrections, the leading EFT corrections as well as the factorised QCD times EFT contributions. For further details see Ref [69]. 5 Constraining dimension-6 EFT coefficients

Monte Carlo
We have implemented the Lagrangian terms of Eq. (3.4) in a HERWIG++ model where the SM matrix elements have been taken from the code HPAIR [10,14,70]. The Monte Carlo event generator is essential to our analysis. This is because a calculation of the total cross section alone cannot account for the change in distributions of momenta and angles that will substantially change the efficiency of the experimental analysis. Thus, even if we do not use their shapes explicitly in extracting constraints in this article, it is essential to have a reliable description of the underlying distributions. One could employ the Monte Carlogenerated distributions to improve the constraints further [71], keeping in mind the fact that operators will mix due to renormalization group running. See for example, Ref. [72]. Throughout this article, we restrict our calculation of the hh production cross section to LO in QCD in the EFT framework (see Fig. 1) and normalize our result to the state-of-theart (NNLO QCD) SM calculation of ∼ 40 fb [17], thus effectively applying a flat K-factor of 2 to all the LO+EFT cross sections. 9 We start by investigating the individual effects on the gluon fusion production cross section, varying one coefficients at a time, while setting all others to zero. The result is shown in Fig. 2, for the LHC running at 14 TeV proton-proton centre-of-mass energy, where we have zoomed-in in the right panel, and shaded the ±10% variation region from the SM value of the cross section. Here and in the remaining article, we employ the MSTW2008nlo_nf4 PDF sets [73]. 10 One can clearly see how deviations from the SM prediction c i = 0 could lead to substantial changes in the total cross section. Unfortunately, the dependence on c 6 is rather mild, whereas the dependence on c t and c g is substantially more pronounced.

Analysis
To accommodate a direct comparison with existing phenomenological analyses, we focus on the process hh → (bb)(τ + τ − ) at the 14 TeV LHC. The specific final state possesses a relatively large branching ratio and manageable backgrounds. This channel has been examined in detail within the SM in [22][23][24]29] and turned out to be particularly promising. 9 We use a rather conservative flat K-factor of 2 as the NNLO calculation of Ref. [17] is not available at present to use for the given parameters that we employ here (i.e. PDF set and scale choices). 10 The cross sections have been verified through an independent implementation directly in HPAIR.
We consider here only the main irreducible backgrounds, arising from tt production with subsequent decays of the W bosons to τ leptons, as well as ZZ and hZ production with (bb)(τ + τ − ) final states, which is sufficient given the other sources of uncertainty. The backgrounds were generated at next-to-leading order in QCD, using the aMC@NLO event generator [74][75][76]. The total cross section for tt was normalised to σ tt = 900 pb [77,78] and the ZZ and hZ NLO cross sections were taken out of the aMC@NLO calculation directly: σ ZZ = 15.25 pb and σ hZ = 0.8329 pb. For realistic description of the final states, parton showering and hadronization were performed using HERWIG++, and the simulation of the underlying event was included via multiple secondary parton interactions [79]. We follow the basic analysis steps as given originally in [22] and as were described in [29]. Here, we assume 70% τ -reconstruction efficiency with negligible fake rate 11 and require two τ -tagged jets with at least p ⊥ > 20 GeV. We require that the di-tau invariant mass, taken from the Monte Carlo truth, reproduce the Higgs mass within a ±25 GeV window, to account for the reconstruction smearing, as done in [22]. We further smear the true di-tau invariant mass by a 20 GeV Gaussian to allow for the contamination of events containing Z → τ + τ − into the di-tau mass window that we consider. We use the Cambridge-Aachen jet algorithm available in the FastJet package [80,81] with a radius parameter R = 1.4 to search for so-called 'fat jets'. We require the existence of one fat jet in the event satisfying the mass-drop criteria as done in the hV study in Ref. [82]. We require the two hardest 'filtered' sub-jets to be b-tagged 12 and to be central (|η| < 2.5) and the filtered fat jet to be in (m h − 25 GeV, m h + 25 GeV), which will remove events with Z → bb decays. The b-tagging efficiency was taken to be 70%, again with negligible fake rate for the sake of simplicity. We require a loose cut on the transverse momentum of the fat jet (after filtering) that satisfies the above criteria, p fat ⊥ > 100 GeV and also apply a transverse momentum cut on the τ + τ − system of equal magnitude, p τ τ ⊥ > 100 GeV. As in [29], we apply additional cuts: ∆R(h, h) > 2.8 and p hh ⊥ < 80 GeV to reject the background even further.
We investigate the effect of the above analysis on events corresponding to different values of the coefficients. 13 On the left panel of Fig. 3 we show the efficiency of the analysis in the case where we set all other coefficients to zero except the labeled one. On the right panel we show the efficiency times the cross section of the EFT point, divided by the SM cross section (σ LO+EFT = 22.3 fb) times the SM point efficiency ( 0.065). This plot can be compared with the actual cross section plot of Fig. 2, where one can observe that the qualitative behaviour of the resulting cross sections does not change, but there is a quantitative change due to the non-uniform effect of the analysis. Note that we have not included the actual branching ratios in Fig. 3, which will have a further effect in determining the significance of a given parameter point.

c 6 -only model
For simplicity, we begin by considering a model which contains only a non-zero coefficient c 6 . In fact, this is the only coefficient that remains unconstrained from data on single Higgs boson production. Setting all other coefficients to zero and varying c 6 corresponds to modifying the SM value of λ, as was done in previous studies. 14 No modifications are expected at the order considered on the Higgs boson decays when including such an operator. We investigate the possible constraints on the c 6 coefficient given the particular model.
Let us assume that in our analysis we obtain S(c 6 ) events for the signal, for a given value of c 6 at a given integrated luminosity. For the background, we obtain B events at the same luminosity. Given that the number of events S and B is large enough, we may assume that they are Gaussian-distributed. The total statistical uncertainty on N (c 6 ) = S(c 6 ) + B is then given by: Therefore, if the relative theoretical uncertainty on the cross section prediction is f th , assuming negligible theoretical uncertainty on the background, 15 an addition in quadrature leads to a total uncertainty of: In particular, in [27] we focused on such a scenario. Moreover, we studied the dependence on the top Yukawa coupling, which (beyond testing consistency with the SM) allows to examine an independent variation of the coefficient of thetLtRh operator in Eq. (3.4), which is possible in extensions of the SM where the Higgs boson is not part of a SU (2)L doublet. 15 The assumption is reasonable since the background higher-order calculations exhibit relatively small variations compared to the hh signal. Moreover, the background prediction can be normalised using a signal-free region. and hence we have that: To obtain the expected constraints, we assume that the underlying theory is indeed the SM, which corresponds to c 6 = 0 in this scenario. In turn, the expected total number of events is N (c 6 = 0). One then needs to compute how many standard deviations δN (c 6 ) away a given N (c 6 ), as predicted from theory, is from N (c 6 = 0). This can be translated into a probability (i.e. a p-value) assuming a Gaussian distribution. The results are presented in Fig. 4 The bounds are weaker for positive c 6 , since this leads to a reduced cross section and thus to a larger statistical uncertainty. The improvement of the 1σ regions is moderate for f th = 0.3 for an increased luminosity, as in that case the uncertainty is dominated by the systematic uncertainty on the theoretical prediction of the signal rates. Improvements on the theoretical description of the process are thus necessary for improving these bounds.

The full model
Generically, one expects several operators to be present. Here we consider the full parameter space, varying the coefficients of Eq. (3.4), as well as c γ and c W , within the currently allowed regions. The parameter c 6 still plays a somewhat distinct role, as it is currently essentially unconstrained, and the information on it will come primarily from multi-Higgs production. Thus, all two-dimensional exclusion planes need to include this parameter. We calculate the p-values in a similar fashion as before: we assume that the standard model is true, i.e. c i = 0 ∀ i, and compare the number of events after the analysis is performed with those expected from the SM, calculating how 'far' they are in terms of the uncertainty δN . We note that, since we are only considering a single observable, the event rate for the given signal point, we do not expect the constraints on the two-dimensional planes to be strong. Nevertheless, the correlations are still interesting and can be improved either by examining several other observables, or optimising the analysis for different points on the parameter space.
To accommodate for the current constraints on the parameters c g , c t , c b , c H , c γ , c W we use the codes HiggsBounds [83] and HiggsSignals [84] on the eHDECAY output. We employ the "effective coupling" mode, where one defines for the decay of the Higgs into the final state X. In particular, the single Higgs cross section is then scaled using the effective coupling to gluons, g hgg . More explicitly, For further details we refer the reader to the HiggsBounds manual [83]. We perform our numerical scan as follows. We take three sets of values for our parameters in each direction of our (7-1)-dimensional parameter grid, covering all coefficients besides c 6 : {0, ±0.01, ±0.1, ±1}, {0, ±0.25, ±0.5, ±0.75} and {0, ±0.15, ±0.2, ±0.4}, while the latter coefficient is allowed to vary in a larger range, i.e. {−2, 3.5}, in steps of 0.5. As mentioned, we assume the coefficients c b and c τ to be equal. In a first step we discard those points where any branching ratio is negative, that is, when the EFT contribution is larger than the SM one, and with opposite sign. 16 The points are later given as input to HiggsBounds, which checks if a point is excluded at the 95% C.L. by collider data (LEP, Tevatron and LHC). This step is numerically fast, and allows to easily discard points where the EFT effects would have generated an excess of events in single Higgs boson studies at the 7 or 8 TeV LHC. The surviving points are fed into HiggsSignals, which performs a multi-dimensional fit to the Higgs observables and outputs a p-value for each point. We discard points which would have given a severe deficit of events in current Higgs data, and thus we keep only points at the 95% C.L., that is, where the p-value corresponds to less than 2 standard deviations from the mean of a Gaussian distribution.
To calculate the allowed two-dimensional regions, we need to marginalize over the remaining dimensions. To accomplish this, for a given point in the (c i , c 6 )-plane, we average over the p-values obtained by varying along the other dimensions within the allowed ranges where N is the total number of points of the parameter space corresponding to the given (c i , c 6 ). To interpret the resulting probability values, we divide by the corresponding one employing the SM values for the coefficients under consideration The 1σ-equivalent contours are thus drawn by finding the iso-curve corresponding top(c i , c 6 ) = exp(−1/2) ≈ 0.607. We do not show the c W coefficient as it turned out to be heavily constrained by current Higgs data, since the branching fraction of the Higgs boson to W W is very sensitive to this coefficient [69]. We first consider the (c H , c 6 )-plane in Fig. 5. The coefficient c H enters all EFT diagrams by changing the Higgs boson wave function in a universal way, and competes with c 6 by reducing the self-coupling contribution in our convention. Evidently, there exists a complicated correlation between c 6 and c H . The coefficient c H could be constrained to lie within |c H | 0.2 or better at the end of the HL-LHC run from hh production, when considering 30% theoretical uncertainty.
We next consider the (c t , c 6 )-plane in Fig. 6. The coefficient c t enters all diagrams that contain top quarks. Points with negative (positive) c 6 and positive (negative) c t are more challenging to exclude, since the coefficients enter in the first line of Eq. (3.15) with the same sign, which leads to a compensation of the effects (see also Fig. 2). It is evident that improving the knowledge on the poorly-constrained 'top Yukawa' c t , entering hh production in various ways, will be helpful to improve the exclusion range for c 6 .
The expected constraints for c g , which adds tree-level couplings of one or two Higgs boson to two gluons, are shown in the (c g , c 6 )-plane in Fig. 7. The constraint on c g is found to be −0.2 c g 0.1 at 3000 fb −1 given that f th = 0.3.
We show the results involving the coefficient c γ in Fig. 8, which enters the process under consideration indirectly, through modification of the branching ratios. The correlation with c 6 is rather modest, with c γ potentially being constrained by this channel at 3000 fb −1 with f th = 0.3 to lie within −0.2 c γ 0.1.
The c b (= c τ ) coefficient is considered in Fig. 9 on the (c b , c 6 )-plane. This coefficient is expected to be sub-dominant in the production due to the assumption of MFV, but it affects the decays of the Higgs boson to bb and ττ , and hence it is relevant to the process we are considering. The coefficient has been already constrained to lie within c b −0.2 by HiggsBounds/HiggsSignals, and thus we show only this region. The coefficient c b can be clearly constrained to the region c b 0.1, even at 600 fb −1 with f th = 0.3.
To obtain a constraint on the c 6 coefficient alone, we marginalize over all the other coefficients. This is done in the same way as prescribed by Eqs. (5.7) and (5.8). The result is shown in Fig. 10, which shows the p-value plotted against c 6 . We also present the 1-sigma contours as black lines.

c 6 − c t model
As a further example, we constrain the only non-zero coefficients to be c 6 and c t , varied in the same regions as in the full model. This can be seen as an approximation towards a scenario where all the other coefficients are constrained substantially by future experiments. As in the previous sub-section, we marginalize over the c t coefficient. The resulting c 6 pvalues are shown in Fig. 11.

Summary of results
We show a summary of constraints on c 6 obtained by all three cases we considered in Table 2.
As explained before, the values take into account the uncertainty due to the other weakly bounded coefficients. As expected, the full model, including the current constraints coming from single Higgs boson measurements, provides the widest 1σ range for the coefficient c 6 , at 3000 fb −1 c 6 (−0.9 c 6 1.7), when f th = 0.3, whereas the c 6 -only model provides, as expected, the most narrow: |c 6 | 0.4 at 1σ at 3000 fb −1 . Given that future measurements  Table 2: A summary of the constraints obtained on the coefficient c 6 at 1σ, at integrated luminosities of 600 fb −1 and 3000 fb −1 at a 14 TeV LHC, assuming that the theoretical uncertainty on the signal rates is 30%.

Conclusions
We have investigated the dimension-6 effective theory description of beyond-the-Standard Model modifications to Higgs boson pair production. Using an implementation within the HERWIG++ Monte Carlo event generator and a realistic analysis, including the description of theoretical uncertainty on the signal rates, we have constructed the possible exclusion regions in a c 6 -only model, a general EFT with all coefficients allowed to vary, and a constrained EFT, marginalising over various coefficients in the latter two cases. We find that at the Large Hadron Collider at 14 TeV, meaningful constraints can be obtained on the parameter space of the EFT, particularly on the c 6 coefficient. These results appear in Table 2. We conclude that, approximately, the hitherto unconstrained c 6 coefficient would be limited to −0.9 c 6 1.7, independently of the values of the other coefficients, and to |c 6 | 1.3 in the simplified c 6 − c t model, at 1σ and for 3000 fb −1 , assuming 30% theoretical uncertainty on the signal rate. This expected bound could be substantially improved for example, by examining other final states originating from hh, improving the experimental analyses by examining differential distributions and future improvements of the theoretical description of the signal process. Without doubt, the process of Higgs boson pair production will be seriously considered as part of the wider programme of constraining the higher-dimensional effective field theory parameter space.

A SM D=6 Lagrangian
In this section we give more details on the Lagrangian terms not explicitly shown in Eq. (2.1). The part of the SM Lagrangian relevant to our study is given by where for simplicity we have only written the Higgs boson couplings to the 3rd generation. The CP-odd effects are given by  In this work we have neglected CP-odd effects for simplicity, but using current Higgs boson data one can already set some constraints on a Higgs CP-odd component (see e.g. [85][86][87] and references therein).
Finally, the four-fermion operators, employing the same basis as used in [58], read where the last two operators are suppressed in scenarios of minimal flavour violation (MFV), and we have restricted ourselves again to the third generation. Furthermore, for simplicity, we omitted operators with leptons. The full set of relevant operators can be found in Table  2 of Ref. [58].

B Electroweak symmetry breaking with dimension-6 operators
Due to the additional terms originating from the dimension-6 operator ∼ |H| 6 , the position of the electroweak minimum changes with respect to the Standard Model prediction. To find the new minimum we consider the potential: which contains the additional interaction. Applying the minimisation condition ∂V /∂|H| 2 = 0, we obtain where we have expanded the square root in the last step. Ignoring the high-energy solution irrelevant to our theory, i.e. taking the + solution, we obtain : Ignoring the Goldstone modes, we can expand the Higgs field |H| in terms of the physical scalar Higgs boson about this minimum, H ∼ (0, (h + v)/ √ 2). We start by examining the terms arising from the SM Lagrangian and the dimension-6 operators that contribute to the kinetic term: where D µ is the covariant derivative that includes all the interactions of the Higgs field with the gauge bosons. After expansion, we arrive at where we have ignored the gauge boson interactions. To canonically normalise the Higgs boson kinetic term and to remove derivative interactions, we consider the following nonlinear transformation: Plugging this into Eq. (B.5), we find the values of the constants a i that cancel all terms except (1/2)∂ µ h∂ µ h: This shift should be performed everywhere in the Lagrangian, and introduces changes in as well as new interactions. We first perform the shift in the terms contributing to the Higgs boson scalar mass. The relevant terms are: where the last term comes from the |H| 6 interaction. We substitute for µ 2 using Eq.
from which we can immediately deduce the Higgs mass: The terms contributing to multi-Higgs production via gluon fusion at the LHC are: L h n = − µ 2 |H| 2 − λ|H| 4 − y tQL H c t R + y bQL Hb R + h.c.
where we have included in the first line the relevant Standard Model terms that will receive corrections from dimension-6 operators. We proceed by deriving the expressions for the triple coupling after expanding about the minimum and canonically normalising via Eq. (B.8). The relevant terms are the same as those that appear in the potential of Eq. (B.1): Expanding this about the electroweak minimum, we get − c 6 λ 8Λ 2 (v 6 + 6v 5 h + 15v 4 h 2 + 20h 3 v 3 + 15v 2 h 4 + 6h 5 v + h 6 ) . (B.14) Omitting terms with h n , n > 4, and constant terms we arrive at It is convenient to also calculate the h 2 , h 3 and h 4 terms as a function of h up to h 4 : Finally, we focus on the fermion-Higgs boson interactions that receive contributions from where f = t, b, ..., with f L,R the left-and right-handed fields, and the first term comes from the SM whereas the second term is a dimension-6 contribution. Substituting in the shift of Eq. (B.8), and keeping terms up to O(Λ −2 ), we obtain The first line gives the expression for the modified fermion mass, and we can re-express Eq. (B.19) in terms of this: The final term that we need to consider is