On the reinterpretation of non-resonant searches for Higgs boson pairs

The detection of production of a pair of Higgs bosons before the end of LHC operation would be clear evidence of New Physics (NP). As searches for non-resonant production of Higgs pairs are being designed it is of particular importance to be able to conveniently present current experimental results in terms of limits in the most 'model-independent' fashion possible. In this article we provide an analytic parametrization of the {\it differential} Higgs-pair production at the LHC in the effective field theory (EFT) extension of the SM. It results from a fit to the theory prediction for the $gg \to hh$ cross section at the 13\,TeV at the LHC. Subsequently the resulting formula is used for a reweighing technique that allows to recast exclusion bounds from ATLAS and CMS HH$\to\gamma\gamma\,b\bar{b}$ searches to any point of the considered EFT parameter space. We demonstrate with a fast simulation of the LHC detectors that with this approach it is possible to cover the continuously of the EFT parameter space, taking correctly into account the efficiencies of signal selections, without the necessity of rerunning a large number of full detector simulations. Finally, the resulting exclusion bounds are confronted to several explicit models such as setups with additional scalars, including 2HDM, vector-like fermions, and minimal composite Higgs models, mapped to the EFT.


Introduction
Examining the production of pairs of Higgs bosons h is a crucial task in the long term strategy of the LHC. Analysing this process offers a unique direct window on the Higgs potential, whose behaviour lies at the heart of many unresolved questions on nature, such as the hierarchy problem, the question of vacuum stability [1][2][3], the feasibility of electroweak baryogenesis [4][5][6], and the potential dynamical trigger for electroweak symmetry breaking (EWSB) [7][8][9] -just to name a few. Of particular importance regarding the search for physics beyond the Standard Model (SM) are differential observables, such as the distribution of the invariant mass of the Higgs pairs m hh (see, e.g., [10][11][12]), since they are especially sensitive to effects of new physics (NP).
As we eventually want to confront corresponding experimental results with theory, we need to assume a theoretical framework, which in our case will be the effective field theory (EFT) extension of the SM, supplementing it with higher dimensional operators. This captures the effects of any yet undiscovered physics beyond the SM, given that it is separated from the latter by a mass gap and in particular is significantly heavier than the scales probed experimentally. In that sense, the setup is 'model independent', providing the most general parametrization of nature, based on what we observe at low energies. On the other hand, any known UV model (residing at scales well beyond the electroweak scale) can be mapped to such an effective Lagrangian, rendering a formulation of experimental results in terms of bounds on its parameters an extremely useful bridge between data and theory. Having at hand analytic formulae which translate between observables in a realistic collider environment and effective couplings will thus be of utmost importance for the experimental results to be properly interpreted.
In this paper we will present on the one hand such a fit of the theory prediction for the Higgspair production cross section in gluon fusion at the LHC in terms of effective couplings, considering both the total cross section as well as the ones differential in m hh and |cos θ * |. At Leading Order (LO) the totality of the HH process is described by those two variables 1 . Beyond that, we will employ the formula to reinterpret exclusion bounds, derived assuming certain benchmark points in parameter space with a given kinematics, for points with modified kinematic distributions. In fact, different distributions in variables, such as m hh , will lead to different efficiencies of selection cuts, which can be accounted for, knowing the distributions in terms of effective couplings. Without such an analytic knowledge, a full monte-carlo (MC) simulation would actually be required for each parameter-space point, which is not feasible for a dense scan of the parameters. Finally, we will use the interpreted results to confront explicit models, mapped to an effective Lagrangian, efficiently with constraints from the LHC. To perform this exercise we collect from literature setups with additional scalars, including two-Higgs-doublet models (2HDM), vector-like fermions, and minimal composite Higgs models and provide a collection of relations between their explicit parameters and effective couplings after EWSB.

Setup
The theory, on which our analysis is based, is the EFT extension of the SM with D = 6 operators O i with coefficients ∼ c i , suppressed by the NP scale Λ 2 (assuming for the moment a linearly realized SM gauge symmetry). 2 The terms relevant for the production of Higgs pairs in gluon fusion to leading order (LO) are contained in the effective Lagrangian [12] Here, L SM is the SM Lagrangian, λ SM = m 2 h /2v is the SM value of the tri-linear Higgs self coupling, with m h ≈ 125 GeV the mass of the Higgs boson, y t is the top-quark Yukawa coupling, and α s is the strong coupling constant. 3 After EWSB (in unitary gauge) H → 1/ √ 2 (0, v + h) T , where v/ √ 2 ≡ |H| ≈ 174 GeV is the vacuum expectation value of the Higgs field, the terms relevant for the gg → hh process become [10][11][12] In fact, the operators in (1) modify the trilinear Higgs self coupling and the top yukawa coupling, parametrized via multiplicative factors κ λ and κ t , where in the SM κ λ = κ t = 1. Beyond that, they induce additional contact interactions between two Higgs bosons and two fermions and between a gluon pair and one or two Higgs bosons, parametrized by c 2 , c g , and c 2g , which vanish in the SM, i.e., c 2 = c g = c 2g = 0. The relations between these effective couplings and the coefficients of the operators in Eq. (1) can be obtained straightforwardly (see, e.g., [12]) and imply the constraint c g = −c 2g . This correlation is however broken in the non-linear realization of EWSB and is not assumed a priori in this work, where we treat all couplings in Eq. (2) as free. In the next section, we will provide a parametrization of the Higgs-pair production cross sections in terms of the effective couplings entering the effective Lagrangian (2).
3 Fit to the (differential) gg → hh cross section We will now describe our fit to the cross section for Higgs pair production in gluon fusion at the 13 TeV LHC in the EFT defined above. 4 We consider the total cross section but as well as the ones differential in m hh and |cos θ * |, with θ * being the polar angle of either of the Higgs bosons with respect to the beam axis (featuring the same |cos θ * | at parton level).
In Ref. [10], the parameters A i , that form a vector A ≡ (A i ) of coefficients of the polynomial Poly(A), were extracted from simulations to describe the total cross section. In this paper, we push the logic further to describe the differential cross section. We slice the kinematic space into bins with sufficient granularity. In each bin, Eq. 3 holds, but the coefficients become bin dependent, A i → A j i , following the matrix-element (ME) integration and the PDF evolution. We define the cross section in a bin j by σ hh,j ≡ σ hh × Frac j , where Frac j is the fraction of events contained in bin j. The differential hh production ratio now becomes We simulate the differential cross section for various points in the five-dimensional parameter space, spanned by {κ λ , κ t , c 2 , c g , c 2g }, to determine A j . At LO accuracy the ME information differential in the invariant mass of the Higgs pair (m hh ) and in |cos θ * | completely defines the 2→2 process. At any subsequent simulation step -after the Higgs decays -the process evolution is independent of the NP at high energies, up to detector level. Therefore the results based on the ME-level information are applicable to reconstruct NP shapes at LO at detector level, provided that the former is known on an event-by-event basis.
After an optimization procedure, we converged to an optimal list of 59 bins in m hh that allowed for a precise reweighing of an ensemble of generated events to any NP shape. Since the distribution in the cos θ * variable is rather flat, only 4 bins are considered 5 . For the binning in m hh we choose 10 GeV intervals up to 700 GeV, where the bulk of the cross section is observed, and use a more coarse binning above. Explicitly, our binning is given by Following the procedure established in [10], the components of A j are extracted by maximizing the likelihood simultaneously for all the coefficients, employing an ensemble of MC simulated samples, scanning the model parameters, i.e., minimizing where the index i runs over different points in the {κ λ , κ t , c 2 , c g , c 2g } parameter space. Here, T (i,j) ≡ σ i hh,M C Frac j i , with σ i hh,M C the total cross section calculated via MC simulation and Frac j i the corresponding bin fraction (also taken from MC), while σ i hh (A j ) ≡ R j hh σ SM hh Frac j SM is the differential cross section parametrization following Eq. 4. Only the statistical (MC) uncertainty on the cross section in each bin, δT (i,j) , is considered. For the SM point we assume a non-zero value of 10 −4 to regularize the likelihood. 5 Here, and in the following, we omit the absolute value signs around cos θ * for brevity.
The simulation of events is performed within the MG5_aMC@NLO framework [15] where the entire event generation process is automated [16]. We make use of UFO model files [17] extracted from the Lagrangian 1, as constructed by the authors of Ref. [18]. For each NP point we simulate 50,000 events, while for the SM benchmark we generate a sample of 13,000,000 events.
Our results for the fit coefficients are compiled in Tables 9-12. To illustrate the bin statistics, we also show the number of events in each bin, both for our total BSM sample (see Section 4) and the SM sample.
In Figure 1 we provide some examples of the fit agreement with the dataset points to selected m hh ranges considering the bin central in cos θ * . Here, we display three parameter space directions, which are κ λ for fixed κ t = 1, κ t for κ λ = 10 −4 and c 2 for κ λ = 1, κ t = 1.5. The errors on the data points are purely statistical, while the coloured lines correspond to the fit result, Eq. (4). It becomes evident from the results that, even if it is true that for every bin one needs only a minimum of N fit = 15 different samples to determine all coefficients, due to limited MC statistics the result of such a procedure would have a limited reliability. Concerning the behaviour of the cross section in different bins, the results confirm the expectation that a change in the scalar self coupling κ λ affects mostly the threshold region, while the (D > 4)tthh contact interaction has a large effect on the high energy tail. Moreover, the pronounced destructive interference present in the cross section in the vicinity of the SM point κ λ = 1 becomes evident via the peaked minimum, visible in the left-most plot.
In Figure 2 we display the value of the coefficients of R j hh for the central cos θ * bin as a function of m hh , grouped in the subspaces given in Table 1. As a first observation, we see that at high invariant masses A 2,5,10 become most important, which coincide with the terms that parametrize genuine higher dimensional effects (containing exclusively c 2,g,2g ), in agreement with expectations. We further note that for the subspace that corresponds solely to a variation of SM parameters (S1,    (4), with the same ratio, as derived directly via MC, for the most central bin in cos θ * (betweem 0 and 0.4) and for four different bins in m hh , which are between 260 GeV and 270 GeV (threshold), 350 GeV and 360 GeV (interference region), 500 GeV and 510 GeV, and 1000 GeV and 1100 GeV (high mass tails). From left to right we show a scan in k λ , k t and c 2 , respectively, while keeping the other parameters at fixed values The error bars have a pure statistical source related with the generated number of events.
with coefficients A 1,3,7 ), the absolute value of the coefficients does not surpass 100, and they are peaked towards low m hh . In particular, the value of A 7 (that controls the triangle-box interference term ∼ κ 2 t κ t κ λ ) is always negative for m hh < 400 GeV. In subspace S2 (coefficients A 2,6,8 ), which features the contributions from the newtthh contact interaction (including interference with SMlike contributions), the coefficients are in general larger compared to those in S1, resulting in stronger variations with the corresponding parameters close to the hh threshold, and show more pronounced m hh tails, as discussed before. Finally, all the remaining terms contain gluon-Higgs contact interactions and as such also feature pronounced contributions at large energies. The Feynman-diagrammatic representation of the different contributions to the cross section can be studied in Figure 1 of [10].
After this general discussion of the behavior of the cross section and the fit coefficients, in the next section we will employ Eq. (4) to reinterpret experimental analyses.  Figure 2: Values of the fit coefficients entering R j hh , differential in m hh . We only display the results for the central cos θ * bin (betweem 0 and 0.4). From left to right, we show the coefficients belonging to subspaces (S1,S2), S4, and (S3, S5, S6), respectively (see Table 1). 6

Recast of ATLAS and CMS measurements
In this section we establish and validate the reweighting method and estimate the present limits on the Higgs couplings in a multi-dimensional space. The constraints are obtained by a recast of preliminary results from the ATLAS and CMS collaborations [19,20] in the HH → γγbb channel, that refer, respectively, to samples of 3.2 fb −1 and 2.7 fb −1 of 13 TeV data. Results with a slightly better sensitivity were published employing 8 TeV data [21][22][23], but we decided to use the preliminary results since they are obtained at the same center-of-mass energy that will be used for the next generation of LHC results. Between the publications of these and the release of this paper, other searches were performed by both collaboration in several channels [21,24,25] and with more luminosity, were the CMS analyses also include results for beyond the SM (BSM) scans in κ λ and κ t .
The choice of the HH → γγbb channel results from the fact that this final state is easy to simulate and to reconstruct using a parametric model of the ATLAS and CMS detectors. Moreover, both analyses were performed with a sequential application of selections (in contrast to a multivariate analysis that is a kind of standard now), making them easy to recast. It is known that this final state is less constrained in regions of the phase-space compared to other final states [21]. Nevertheless, since our goal is the proof of principle of the recast technique, we are not hunting for the latest and best limit. In the case of more sophisticated (multivariate) analyses, a reinterpretation of the results is more subtle, however still possible if samples exist for all relevant kinematic configurations, as provided in [10] (see discussion at the end of Section 5).

Signal simulation and reweighing
We employ Eq. 4 to reweight a set of base events to different points in the parameter space. As base events, we use the 12 benchmarks defined in [10,26,27] and for each of them we simulate 100,000 events. Due to the construction of the benchmarks, the resulting ensemble S BSM of 1,200,000 events populates well the EFT space we are probing. In addition 100,000 SM-like events are generated for comparison but kept out of S BSM .
The hard scattering events are generated at LO accuracy using the model described in section 3. Following Ref. [26] we normalize the events to a total cross section with a k-factor is applied to include corrections up to NNLO (augmented with NNLL re-summation) in QCD (σ SM hh = 33.5 ± 1.4 fb [27]). Note that the normalization only becomes relevant to predict the absolute number of events, as used in Section 5. In fact, up to NNLO, the k-factors are found to be mostly flat in kinematic variables as well as as a function of the EFT parameters (in the bulk of the parmeterspace) [28] and this suggests that the distributions are well described by performing the reweighing with LO information and employing flat k-factors to correct for QCD corrections at the end is a reasonable approximation. Furthermore, we used the NLO set of the PDF4LHC parton densities [29-32] as well as the parton shower and hadronization infrastructure of the Pythia 8 package [33]. Finally, we simulate the response of an LHC-like detector by using Delphes 3 [34].
To reweight simulated events from S BSM to any other point in parameter space (such as the SM point, which we will use for validation), we apply an event-by-event weight where m i hh and cos θ * i refer to the ME level variables corresponding to event i, and N i is the number of events in bin j that contains m i hh and cos θ *  lines, R hh (m i HH , cos θ * i ) is given by R j hh (Eq. (4)) for this bin. It is our choice to normalize the weights in such a way that the resulting BSM distribution is normalized to one if no cut is applied and to the signal efficiency if any selection is applied. In this manner we decouple the shape information from the total cross section information.
The normalization coefficient C norm is obtained as and is used to unitarize the reweighed shape. In the ideal case C norm = 1. Nevertheless, fluctuations in the samples used to derive the weights can induce a departure from unitarity. For the SM shape C norm = 1, since the bin-by-bin fits are performed fixing the SM point and statistical errors occur only for BSM points where for some samples there are bins with inevitably low statistics (for example in the high m HH bins of threshold-like EFT points). After the renormalization procedure we observe a good agreement within the statistical uncertainties between MC simulation and reweighed samples, as can be seen in the next section.

The analyses
We consider the number of events N obs observed by each experiment in a signal region and B the number of the expected background events. In the ATLAS documents both values are directly provided in the text, while in the CMS document they are extracted from figures imposing a M γγ window identical to the ATLAS one. Both experiments found that the observed events are statistically compatibles with the expected SM backgrounds. Therefore, we can use this observation to derive exclusion limits on the maximal allowed number of (expected) events in a NP scenario N BSM that could be compatible at 95% CL with the observed data. We build a Poissonian statistical likelihood where S is the number of Signal events. Bayesian hypothesis with a uniform prior on S allows to invert the likelihood to obtain a posterior as a function of S, normalized to unity between S = [0, ∞]. The likelihoods are shown in Figure 3. We integrate then the likelihood from 0 to 95% to obtain the 2 standard deviation confidence level value for N BSM . The number of observed, expected and excluded events are provided in Table 2. In this table we also include the number of events that would be excluded at L=100 fb −1 , projected assuming N obs (100 fb −1 ) = B(100 fb −1 ) = L * B/L , where L is the data luminosity currently used.
To translate the limit on N BSM into the limit on the BSM cross section σ BSM we need to know the selection efficiency within the signal box and the luminosity L: σ BSM = N BSM /(L ). For both analyses we can extract (or estimate) SM , i.e., the efficiency for the SM-like hh production, from the information provided in the respective papers. What is not known a priori is the efficiency 8  for a given BSM model with the associated choice of parameters. To estimate this, we employ a mock-up of the ATLAS and CMS analyses with Delphes 3 . The resulting efficiency D is then rescaled (if needed) by a fudge factor f SM = SM / SM,D to get an estimate of the real efficiency This factor takes into account the limited information we benefit from to emulate the analysis and the approximations made in Delphes 3 to describe the detector response, such as the efficiency of detecting photons or tagging b-jets. Note that this step could be avoided if the experimental collaborations would make public the full information (including the generation level m hh and |cos θ * | variables) on the events that pass the signal selection.
In Table 3 we summarize the set of selections of both ATLAS and CMS searches. It is important to notice that the same Delphes 3 card was used for both experiments. This approximation is valid for photons in the energy range considered in hh production since both experiments are more than 90% efficient to tag them. The difference in fiducial acceptance is emulated at the level of selections. For jets, both experiments use the anti-k T algorithm [35] with radius parameter 0.4. Their acceptance and calibrations are similar at those energies. The major difference comes from the b-tagging working point (WP). Since both experiments require both jets to be b-tagged we just rescale the typical efficiency of the WP used in Delphes 3 within our p T and η range (63% -see Ref. [36]) to the efficiency of the WP declared in the public notes, i.e., 85% for ATLAS and 78% for CMS. Finally, for the projections to 100 fb −1 given below, we assume the observed number of events to match the expected and no improvements in the analysis, which is a very conservative assumption.
After the reweighting according to the b-tagging properties and the selections of Table 3, we find the efficiencies shown in Table 4. It appears that for the ATLAS analysis our efficiencies have to be rescaled by f SM = 1.41. Beside the usual suspects described above to explain the size of the factor, we suspect that the photon efficiency is slightly better in the real ATLAS analysis than in the Delphes 3 card, but the exact amount is hard to guess from the details given in the conference note.
In the CMS case the efficiency of the 2 b-tag and 1 b-tag categories are only provided together and no explicit information on how it is distributed between them is given. One may notice that a similar categorization was used at 8 TeV where the events appeared to be evenly split between the two categories, incidentally [23]. We apply the same assumption to the 13 TeV analysis and split the efficiencies evenly between the two categories ( SM,CMS ≈ 10% per category).

Validation of the analytical reweighing
In Figure 4 we show the comparison of the actual simulation and the reweighting for an SM-like signal for 3 main reconstructed kinematic variables describing the hh system: m X ≡ m γγbb − m bb − m γγ + 250 GeV, cos θ * γγbb and p T,γγbb . The first variable was shown to be the best estimate of m hh in Ref. [20]. The statistical errors per bin are shown only for the plain simulation, since the shapes derived by reweighting are constructed from 12 times as much bare events and thus expected to have less statistical fluctuations. In general we find a very good agreement between the distributions. Comparisons for different benchmarks provide similar picture and are given in Appendix B. Moreover, in Figure 5 we show reweighted distributions in m X , after ATLAS-like selections, for BSM points not within the set of (simulated) benchmarks. Beyond points featuring non-vanishing contact interactions (cg=-c2g=1 and c2=1), we also consider the maximum boxtriangle interference (kl=2.4), thereby exploring three different kinematic regimes. In all the cases we see smooth shapes.
Finally, the efficiencies of different benchmark points obtained from the MC simulation are given in Table 5 ( M C ), together with the reweighted efficiencies as obtained after applying the renormalization via C norm . The agreement of both efficiencies is remarkable, specially in the points where EFT is linear (where the parameter space scan used for the fit is more dense, see table 1).

Model dependent interpretations
We illustrate the method described in Section 4 reinterpreting EFT bounds, applying them to concrete NP setups. We consider the explicit models collected in Table 6.
In Table 7 we provide the relations between the effective couplings for the various models   under consideration, arising after integrating out the NP at tree level, where we express dependent couplings in terms of independent ones, the latter given in the second column. The relation with the fundamental parameters of the models is given in Table 8, which contains equivalent information.
Note that the fact that we integrate out NP at tree level might appear problematic, given that the effective couplings in Table 7 enter Higgs-pair production only at one-loop and, in principle, potential loop induced contact interactions with gluons (c g,2g ) could enter the process at tree-level, thus lifting a loop-suppression in integrating out NP. However, since the scalars considered in models 1-7 are assumed to be color neutral, no such effective interactions are in fact induced at the one-loop level. The same is true for integrating out vector-like leptons (model 9). The vector-like quarks (model 8), on the other hand, do not couple to the Higgs boson without involving further SM fermions, thus also not generating Higgs-gluon contact interactions at one-loop. Table 7: Correlations between effective couplings in various explicit models. While the total number of free parameters is invariant, when there was freedom which parameters to treat as free, we always chose κ t , expressing the other couplings in terms of the latter.

Model Free Parameters
The situation becomes more subtle for the composite Higgs setups (models 10 and 11). Here, the Model Fund. Parameters Table 8: Effective couplings in terms of physical parameters of various explicit models. Here, α is the mixing angle between the two scalars, while β = arccos(v 1 /v) is the arccosine of the ratio of the vev of the (first) doublet and the electroweak vev v ≈ 246 GeV, and we defined s x ≡ sin x (t x ≡ tan x). A common mass m H is assumed for the heavy scalars, besides for H + and A in models 3-6, with masses m H + , m A . Beyond that, m 2 is the coefficient of the triple-singlet coupling and λ α that of the bi-quadratic scalar term, while Z 6 multiplies |H 1 | 2 H † 1 H 2 in the 2HDM. Moreover, M T and M E are the masses of the heavy vector-like quark and lepton, respectively, and λ T t , λ E are the coefficients of their (Yukawa-type) couplings with the SM fermions, mediated by the Higgs. Finally, ξ ≡ v 2 /f 2 parametrizes the composite Higgs non-linearity, with f the Pseudo-Goldstone decay constant. See references given in Table 6 for more details.
potential effect of loop-induced Higgs-gluon contact interactions due to integrating out fermionic resonances appearing in the models (that could become relevant in parts of the parameter space) cancels with additional corrections to the Yukawa couplings generated by the very same same fields. Thus, considering only the Higgs-non-linearities that lead to the anomalous couplings as given in Table 7 (see also Table 8) leads effectively to a correct description. In summary, the effective couplings given in Table 7 provide an appropriate description of all models at hand to leading approximation, given that the NP is heavy such that the EFT framework is valid (see, e.g., [13]).
Using the map from the model parameters to the anomalous couplings in Tables 7 and 8, we determine the sensitivity of the LHC analyses to the independent model parameters (which differs between the various models due to the different correlations). The reweighted differential information is used to determine the expected number of events that populate the signal region in each analysis and a parameter point is excluded if it exceeds the upper limits derived in Section 4. Given the low amount of data in the analyses considered, together with the inherent difficulty of the channel, we do not expect to provide competitive limits on the models at the current stage, but rather present these results as an academic exercise to demonstrate the application of our tool in the future.  We begin with the limits obtained for models with additional scalars in the κ λ − κ t plane, given in Figure 6. In the left plot, the cases of a real scalar singlet with explicit Z 2 breaking, a complex triplet scalar, and the 2HDM are compared with the scenario where only these two effective couplings are allowed to vary and the others are set to zero 8 . Note that in each of the models considered, only κ t between 0 and 1 are permitted 9 . In these examples, the presence of c 2 (for κ t = 1) leads to mild variations in the final sensitivity in this plane. The impact of c 2 on the bounds is further quantified in the right plot, where we show the effect of setting this coefficient to ±1, which shifts the contours left or right -in agreement with the tendencies observable in the left plot.
Moving to the models which can be described with one parameter (combination), we show in Figure 7 the fiducial cross section as a function of the free (fundamental) parameter for the MCHM (left), real singlet scalar with spontaneous Z 2 breaking (center), and singlet vector-like fermion (right) models. For comparison, we also provide current exclusion limits from ATLAS and CMS as well as projections to 100 fb −1 , that will exclude the part of the parameter space corresponding to a larger fiducial cross section.
For both MCHM scenarios, the signal efficiency is observed to be rather flat as the kinematics are close to the SM-like case. The sensitivity thus is mostly determined via a simple rescaling of the total rate. The current analyses are sensitive to values of ξ ∼ 0.9 in the case of the MCHM 5 and values of ξ ∼ 0.2 can be probed with 100 fb −1 , while no bound is obtained for the MCHM 4 . We further operators generated in the models at hand that can modify the Higgs boson branching ratios (in particular h → bb). As the Higgs boson is a narrow particle those effects do not change the hh kinematics and the effect is secondary with respect with the scope of this paper. 8 That has been considered by the CMS collaboration 9 While a large depletion seems only viable in the presence of cancellations with other NP contributions to single Higgs boson production and decay, here we focus mostly on hh production and only take into account rigorous theoretical constraints at this stage, leaving a combined phenomenological analysis for future work. observe a similar pattern in the case of the singlet, where the kinematics are not greatly affected and no bound is obtained for mixing angles cos α 0.5.
In the case of vector-like fermions, the free parameter scales with the Yukawa mixing over the mass of the new states and a decoupling is exhibited with the latter approaching infinity, as expected. As shown in Table 7, the vector-like quark model does not modify κ λ . Combining this with the sensitivity information of Figure 6, it is not surprising that (expected) limits are rather weak, stemming mostly from a non-vanishing c 2 and residing in the region of large couplings (λ T t 1) and/or small masses (M T 1 TeV). The vector-like leptons, on the other hand, do modify κ λ and also identically affect κ t . Since c 2 is unchanged, they correspond exactly to the benchmark scenario considered in the experimental analyses (with the additional constraint κ t ≡ κ λ ) and the limits are also rather weak.
Finally, Figure 8 displays the limits obtained for the models with more than one free parameter, this time in terms of the different 'fundamental' model parameters (and including the real scalar triplet): the real singlet with explicit Z 2 symmetry breaking (left), the real/complex triplet (center) and the 2HDM (right). Where possible, we reduce the parameter space to two degrees of freedom: In the triplet models we set the heavy scalar masses to be equal while for the real singlet we can define an effective trilinear coupling In the 2HDM, it is not trivial to reduce the (t β , Z 6 , m H ) parameter space. However, Z 6 is related to the usual alignment parameter, cos β−α and the neutral scalar masses by [39] v which we use to fix a relation between the latter, showing predictions for two fixed values of Z 6 = ±2.5. 10 The current constraints on the real singlet lie in regions of large λ eff of order 5-10 and of sizable mixing with the Higgs. The bounds will however improve significantly with increasing luminosity. For the triplets, the analyses are sensitive to large mixings for masses of a few hundred GeV while increasing the heavy scalar mass leads to an increased sensitivity, including also smaller mixing angles β. This subtle behaviour is due to the fact that m 2 H is proportional to the size of scalar quartic interactions (times the electroweak vev squared), such that the coupling strength increases with the scalar mass. We thus cut off the parameter space at m 2 H ∼ 4πv 2 ∼ 1 TeV, in order to remain at a reliable (perturbative) behaviour of the theory. Finally, in the case of the 2HDM, a sensitivity to scalar masses in the TeV region is only possible for small tan β, which amplifies the corrections to κ t and c 2 .

Mapping between smooth scan and shape benchmarks
Most of the times the experimental searches are not easy to recast as soon as shape analysis or multivariative analysis is used. In this case we may ask ourselves if at some level of confidence we can use results obtained for shape benchmarks for a general study of the parameter space. In this chapter we test this premise on the scans presented in the last chapter. We verify how the same criteria used to define the shape benchmarks as representatives of clusters of similar shapes from  Table 6 and 95% CL exclusion bounds. See text for details.  Table 6, as derived in the EFT framework.
a large parameter space scan can help to predict the closest experimental limit (provided there is a list of results for the shape benchmarks), and how this limit can approximate the "real" limit . As a quick reminder, in Ref. [10] a two sample Test Statistic (TS) was defined to order the degree of similarity between two samples. The log-likelihood ratio function of the hypothetical case in which the two samples under test share the same parent distribution is the product over the bins of the probability to observe n i,1 and n i,2 event counts in bin i from the two samples S 1 and S 2 and can be written as: This quantity is constructed in a manner that it is "χ 2 distributed" [44,45] and therefore can be directly used as an ordering parameter to decide between pairs of test samples which of them are the most likely to be compatible with the same parental distribution. In other words, the values T S ij and T S kl obtained respectively by testing the compatibility of samples ij and kl are suitable to determine if samples S i and S j are more similar to each other than are samples S k and S l : this is the case if T S ij > T S kl . 16 One solution for constructing analyses optimal in continuous scans is, for example, the use the T S as quantity to decide, given a test sample to which shape benchmark this one is most similar to and check if this prediction corresponds to the closest experimental limit. At this point we should emphasise that we are after the closest experimental limit. The real experimental limit, that would be directly derived for the test point, would be a bit different from the one from the closest shape benchmark.
Since we are considering cut-and-count examples the "experimental limits" can be well identified with the signal efficiency in the signal region. As the CMS-like analysis contains a cut in the m γγbb variable the signal efficiency is more sensitive to the BSM physics, therefore we will only do this mapping exercise for the CMS-like case. In figures 9 and 10 we show the efficiency maps for the cases of the two 2HDM scans and the triplet and singlet extension cases. Of course the quality of color interpolation on the figures depend on the density of points inspected. To each inspected point a marker is superimposed, symbolizing the closest shape benchmark according to eq. 12. The first fact to notice is that as suggested in [10], the regions belonging to the same benchmark tend to enclose fully connected regions of parameter space, that we will call islands. We also notice that within those fully connected regions the closest benchmark indeed has a signal efficiency that is the closest to the true one. This conclusion is not so precise near the boundaries between islands. This approach has multiple applications: If used by the theory side it is an approximation of reality, that can be employed to obtain a first estimate of the constraints on a specific model without the need to expensive MC simulations or recast. This usage would hold for the case of including further BSM effects such as the inclusion of additional EFT operators (for example, the chromomagnetic term [46]) or the departure from the EFT framework towards explicit inclusion of particles in the loops (for example [47,48]). On the experimental side usage, we can imagine a situation where it is advantageous for the discovery potential of a hypothetical BSM signal to construct different selections (or multivariate variables) for the different shape benchmarks. In this case it would be unclear how to extend the interpretation of the search to a smooth scan to decide which set of cuts for use to a point that is not a shape benchmark, and the T S-test offers a solution for keeping the best search sensitivity also for the case of smooth parameter space scans.

Implementation in Rosetta
Although the numerical values of the coefficients that can be used to construct the signal predictions are given in appendix A, for convenience, the input and algorithm to perform sample reweighing are implemented as a module of the Rosetta package [49]. Rosetta provides a framework for EFT basis translations along with a suite of modules ranging from the calculation of Higgs branching fractions via an interface to eHDECAY [50] to the compatibility of a given parameter space point to electroweak precision data and Higgs signal strength measurements. All modules can be exploited in a basis-independent way once a basis definition has be implemented in the package. The package may be found on http://rosetta.hepforge.org/.
The results of this paper are implemented as in extension of the existing dihiggs module [51]. In this module we had implemented the example of usage of the elements of the machinery that do not require expensive Monte Carlo simulations, namely the mapping of parameter space points into shape benchmarks done in the last section. From this implementation it is straightforward for the user to extend it to the case where the calculation of event weights is of interest. After downloading the package, one can invoke the dihiggs interface by calling the command line executable,

>> bin/rosetta dihiggs [OPTIONS]
The package is designed to receive SLHA formatted parameter cards specifying the coefficients of an existing basis implementation. The required translations are then performed to obtain the anomalous Higgs couplings parameters relevant for pair production and decay. Alternatively, one can directly specify the values of the anomalous coupling parameters via the self-explanatory -kl, -kt, -c2, -cg and -c2g optional arguments. The input parameter card has therefore been made optional, to be specified by the -param_card option.
Along with computing the inclusive Higgs pair production cross section, the module will also compute the correct Higgs branching fractions for that point in parameter space using either eHDECAY or an internal interpolating function. Using the results of this paper, we have also added the functionality to return the closest benchmark according to the test defined in formula 12. Finally, the computation of the higgs-pair production process is promoted to differential level, via an internal function. Inside interfaces/dihiggs/AnalyticalReweighter.py the user can find a function weight(variables), where a matrix of weights is calculated from vectors containing the generation level variables m HH and cos θ * according to formula 7. Support for implementing the same algorithm in C++ fashion may be given by directly emailing us. 18

Conclusions
Understanding the properties of the Higgs potential is of utmost importance to eventually answer one of the key questions in particle physics, namely what is the origin of EWSB. Examining the production of Higgs pairs is a crucial experiment to achieve this goal, with a focus both on the total cross section, as well as on distributions, which contain valuable information on the nature of potential NP. In this article we presented semi-analytic results for the distributions of the Higgs-pair production cross section in a well defined and maximally general parametrization of nature, as it appears at low energies, employing EFT. This allows us to express potential deviations from the SM in a consistent way as coefficients of effective operators and constraints on these operators will provide us a guidance on how nature could look at shortest distances.
Furthermore, we employed the formula to recast exclusion bounds, derived assuming certain benchmark points in parameter space with a given kinematics, to points with modified kinematic distributions. The presented method is crucial to cover the full EFT parameter space, taking correctly into account the efficiencies of signal selections.
Finally, the results presented are also useful to confront explicit models, mapped to an EFT, efficiently with constraints from the LHC, providing a bridge between (explicit) theories and data. We demonstrated this procedure, using recent ATLAS and CMS results, for various NP setups, like models with additional scalars, including 2HDM, vector-like fermions, and minimal composite Higgs models, delivering also a dictionary between their explicit parameters and effective couplings after electroweak symmetry breaking.

19
A Numeric tables with the coefficients Tables 9 to 12 contains the numeric values of the formula 3. The first columns of the tables stand for the number of the bin, followed by the mean m hh and cos θ * of that bin. It is also displayed the number of events found in that bin for a 13,000,000 events SM sample and for the 1,200,000 events of the reweighting sample S BSM used in this paper. The same informations can be found visually in figures 11 (for the number of events/bin) and 12 (for the values of the coefficients on the not central cos θ * bins).  Figure 12: Values of the fit coefficients entering R j HH , differential in m HH . We only display the results for the central cos θ * bin (betweem 0 and 0.4). From left to right we show respectivelly the coefficients relative to subspace S1 and S2, S4 and finally S3, S5 and S6 (see table 1