Signs of Composite Higgs Pair Production at Next-to-Leading Order

In composite Higgs models the Higgs boson arises as a pseudo-Goldstone boson from a strongly-interacting sector. Fermion mass generation is possible through partial compositeness accompanied by the appearance of new heavy fermionic resonances. The Higgs couplings to the Standard Model (SM) particles and between the Higgs bosons themselves are modified with respect to the SM. Higgs pair production is sensitive to the trilinear Higgs self-coupling but also to anomalous couplings like the novel 2-Higgs-2-fermion coupling emerging in composite Higgs models. The QCD corrections to SM Higgs boson pair production are known to be large. In this paper we compute, in the limit of heavy loop particle masses, the next-to-leading order (NLO) QCD corrections to Higgs pair production in composite Higgs models without and with new heavy fermions. The relative QCD corrections are found to be almost insensitive both to the compositeness of the Higgs boson and to the details of the heavy fermion spectrum, since the leading order cross section dominantly factorizes. With the obtained results we investigate the question if, taking into account Higgs coupling constraints, new physics could first be seen in Higgs pair production. We find this to be the case in the high-luminosity option of the LHC for composite Higgs models with heavy fermions. We also investigate the invariant mass distributions at NLO QCD. While they are sensitive to the Higgs non-linearities and hence anomalous couplings, the influence of the heavy fermions is much less pronounced.


Introduction
The LHC Higgs data of Run 1 suggest that the scalar particle observed by the LHC experiments ATLAS and CMS in 2012 [1,2] is compatible with the Higgs boson of the Standard Model (SM). The non-vanishing vacuum expectation value (VEV) v of the SU (2) Higgs doublet field φ in the ground state is crucial for the mechanism of electroweak symmetry breaking (EWSB) [3]. It it is induced by the Higgs potential Introducing the Higgs field in the unitary gauge, φ = (0, In the SM the trilinear and quartic Higgs self-couplings are uniquely determined in terms of the Higgs boson mass M H = √ 2λv, with v ≈ 246 GeV. The experimental verification of the form of the Higgs potential through the measurement of the Higgs self-couplings is the final step in the program aimed to test the mechanism of EWSB. The Higgs self-couplings are accessible in multi-Higgs production processes [4][5][6][7]. While previous studies [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] showed that the probe of the trilinear Higgs self-coupling in Higgs pair production should be possible at the high-luminosity LHC, although it is experimentally very challenging, the quartic Higgs self-interaction is out of reach. The cross section of triple Higgs production giving access to this coupling suffers from too low signal rates fighting against a large background [5,7,27]. The relations in Eq. (3) do not hold in models beyond the SM (BSM), and this would manifest itself in the Higgs pair production process. In general, however, new physics (NP) not only affects the value of the Higgs self-coupling, but also other couplings involved in the Higgs pair production process. 1 An approach that allows to smoothly depart from the SM in a consistent and model-independent way is offered by the effective field theory (EFT) framework based on higher dimensional operators which are added to the SM Lagrangian with coefficients that are suppressed by the typical scale Λ where NP becomes relevant [29][30][31][32][33]. These higher dimensional operators modify the couplings involved in Higgs pair production, such as the trilinear Higgs self-coupling and the Higgs Yukawa couplings. Additionally they give rise to novel couplings, like a 2-Higgs-2-fermion coupling, that can have a significant effect on the process. While the trilinear Higgs self-coupling has not been delimited experimentally yet 2 , the Higgs couplings to the SM particles have been constrained by the LHC data and in particular the Higgs couplings to the massive gauge bosons. An interesting question to ask is, while taking into account the information on the Higgs properties gathered at the LHC, if it could be that despite the Higgs boson behaving SM-like, we see NP emerging in Higgs pair production? And if so, could it even be, that we see NP before having any other direct hints e.g. from new resonances or indirect hints from e.g. Higgs coupling measurements?
Previous works have applied the EFT approach to investigate BSM effects in Higgs pair production. 3 A study of the effects of genuine dimension-six operators in Higgs pair production can be found in Ref. [65]. Anomalous couplings in Higgs pair production have been investigated in [66][67][68][69]. In [70][71][72] the EFT was applied to investigate the prospects of probing the trilinear Higgs selfcoupling at the LHC. Reference [73] on the other hand addressed the question on the range of validity of the EFT approach for Higgs pair production by using the universal extra dimension model.
The dominant Higgs pair production process at the LHC is gluon fusion, gg → HH, which is mediated by loops of heavy fermions. It can be modified due to NP via deviations in the trilinear Higgs self-coupling, in the Higgs to fermion couplings, via new couplings such as a direct coupling of two fermions to two Higgs bosons, new particles like e.g. heavy quark partners in the loop, or additional (virtual) Higgs bosons, splitting into two lighter final state Higgs bosons. The purpose of this paper is to address the question of whether it will be possible to see deviations from the SM for the first time in non-resonant Higgs pair production processes by considering explicit models. It has been found that large deviations from SM Higgs pair production can arise in composite Higgs models, which is mainly due to the novel 2-Higgs-2-fermion coupling [36,74]. In this paper, we will hence focus on this class of models. We assume that no deviations with respect to the SM are seen in any of the LHC Higgs coupling analyses, i.e. that the deviations in the standard Higgs couplings due to NP are below the expected experimental sensitivity, for the case of the LHC highenergy Run 2 and for the high-luminosity option of the LHC. Additionally, we assume that no NP will be observed in direct searches or indirect measurements. The prospects of NP emerging from composite Higgs models for the first time in non-resonant Higgs pair production from gluon fusion are analyzed under these conditions. Our analysis is complementary to previous works [46,75], which focused on deviations in Higgs pair production due to modifications in the trilinear Higgs coupling. In Ref. [75] the question is investigated on how well the trilinear Higgs coupling needs to be measured in various scenarios to be able to probe NP. The main focus of Ref. [46] is on how to combine a deviation in the trilinear Higgs coupling with other Higgs coupling measurements to support certain BSM extensions.
Gluon fusion into Higgs pairs exhibits large QCD corrections. In Ref. [4], the next-to-leading order (NLO) QCD corrections were computed in the large top mass approximation and found to be of O(90%) at √ s = 14 TeV for a Higgs boson mass of 125 GeV. The effects of finite top quark masses have been analyzed in [37,[76][77][78][79][80]. While the m t → ∞ approximation exhibits uncertainties of order 20% on the leading order (LO) cross section at √ s = 14 TeV for a light Higgs boson [74,81,82] and badly fails to reproduce the differential distributions [9], the uncertainty on the K-factor, i.e. the ratio between the loop-corrected and the LO cross section, is much smaller due to the fact that in the dominant soft and collinear contributions the full LO cross section can be factored out. The next-to-next-to-leading order (NNLO) corrections have been provided by [83][84][85] in the heavy top mass limit. The finite top mass effects have been estimated to be of about 10% at NLO and ∼ 5% at NNLO [86]. Soft gluon resummation at next-to-leading logarithmic order has been performed in [87] and has been extended recently to the next-to-next-to-leading logarithmic level in [88]. First results towards a fully differential NLO calculation have been provided in [78,80]. For a precise determination of the accessibility of BSM effects in gluon fusion to a Higgs pair, the NLO QCD corrections are essential and need to be computed in the context of these models. They have been provided in the large loop particle mass limit for the singlet-extended SM [58], for the 2-Higgs-doublet model [48] and for the MSSM [4,63]. 4 In the same limit, the NLO QCD corrections including dimension-6 operators have been computed in [89]. In this work, we calculate for the first time the NLO QCD corrections in the large loop particle mass limit for models with vector-like fermions such as composite Higgs models.
The paper is organized as follows. In Section 2 we briefly introduce composite Higgs models. In section 3 we present the NLO QCD corrections to the gluon fusion process in the framework of composite Higgs models including vector-like fermions. In the subsequent sections we analyze whether a possible deviation from the SM signal could be seen or not at the LHC Run 2 with an integrated luminosity of 300 fb −1 and/or the high-luminosity LHC with an integrated luminosity of 3000 fb −1 for different models: in section 4 for the composite Higgs models MCHM4 and MCHM5, and in section 5 for a composite Higgs model with one multiplet of fermionic resonances below the cut-off. In section 6 we discuss the invariant mass distributions with and without the inclusion of the new fermions. We conclude in section 7.

Composite Higgs Models
In composite Higgs models the Higgs boson arises as a pseudo-Nambu Goldstone boson of a strongly interacting sector [90][91][92][93][94][95][96]. A global symmetry is broken at the scale f to a subgroup containing at least the SM gauge group. The new strongly-interacting sector can be characterized by a mass scale m ρ and a coupling g ρ , with f = m ρ /g ρ . An effective low-energy description of such models is provided by the Strongly Interacting Light Higgs (SILH) Lagrangian [97], which, in addition to the SM Lagrangian, contains higher dimensional operators including the SM Higgs doublet φ to account for the composite nature of the Higgs boson. Listing only the operators relevant for Higgs pair production by gluon fusion, the SILH Lagrangian reads 5 with the Yukawa couplings y q = √ 2m q /v (q = u, d), where m q denotes the quark mass, λ the quartic Higgs coupling and α s = g 2 s /(4π) the strong coupling constant in terms of the SU (3) c gauge coupling g s . 6 Here Q L denotes the left-handed quark doublet. The effective Lagrangian accounts for several effects that can occur in Higgs pair production via gluon fusion in composite Higgs models: a shift in the trilinear Higgs self-coupling and in the Higgs couplings to the fermions, a novel coupling of two fermions to two Higgs bosons and additional new fermions in the loops. The latter effect is encoded in the effective operator with the gluon field strength tensors G µν coupling directly to the Higgs doublet φ. While the SILH Lagrangian Eq. (4) is a valid description for small values of ξ = (v/f ) 2 , larger values require a resummation of the series in ξ. This is provided by 4 Reference [63] also shows how the provided results can be adapted to the Next-to-Minimal Supersymmetric extension of the SM. 5 We have not included the chromomagnetic dipole moment operator which modifies the interactions between the gluons, the top quark and the Higgs boson and can be expected to be of moderate size [98]. 6 The relation between the coefficients c and the coefficients c in Eq. (2.1) of Ref. [89] is cx = cxξ (x = H, 6), cu = cuξ and cg = α2/(16π)y 2 t /g 2 ρ cgξ with ξ = v 2 /f 2 and α2 = √ 2GF m 2 W /π in terms of the Fermi constant GF and the W boson mass mW .

MCHM4
MCHM5   Table 1 we report the modifications of the Higgs couplings to the SM particles with respect to the corresponding SM couplings in the SILH set-up and in the MCHM4 and MCHM5. The last two lines list the novel couplings not present in the SM, i.e. the 2-Higgs-2-fermion coupling and the effective single and double Higgs couplings to a gluon pair, as defined in the Feynman rules derived from the SILH Lagrangian, where k 1,2 denote the incoming momenta of the two gluons g a µ (k 1 ) and g b ν (k 2 ). The effective gluon couplings are not present in MCHM4 and MCHM5.
In composite Higgs models fermion mass generation can be achieved by the principle of partial compositeness [101]. The SM fermions are elementary particles that couple linearly to heavy states of the strong sector with equal quantum numbers under the SM gauge group. In particular the top quark can be largely composite. But also the bottom quark can have a sizeable coupling to heavy bottom partners. For gluon fusion this not only means that new bottom and top partners are running in the loops but mixing effects also induce further changes in the top-and bottom-Higgs Yukawa couplings. In addition to the MCHM4 and 5 models involving only the pure non-linearities of the Higgs boson in the Higgs couplings, we consider a model with heavy top and bottom partners based on the minimal SO(5) × U (1) X /SO(4) × U (1) X symmetry breaking pattern. The additional U (1) X is introduced to guarantee the correct fermion charges. The new fermions transform in the antisymmetric representation 10 of SO(5) in this model MCHM10, given by with the electric charge-2/3 fermions u, u 1 , t 4 and T 4 , the fermions d, d 1 and d 4 with charge −1/3, and the χ, χ 1 and χ 4 with charge 5/3. The coset SO(5)/SO(4) leads to four Goldstone bosons, among which three provide the longitudinal modes of the massive vector bosons W ± and Z, and the remaining one is the Higgs boson. The four Goldstone bosons can be parameterized in terms of the field with the generators Tâ (â = 1, ..., 4) of the coset SO(5)/SO(4) The generators of the SU (2) L,R in the fundamental representation read (a, b, c = 1, 2, 3, i, j = 1, ..., 5), The non-linear σ-model describing the effective low-energy physics of the strong sector is given by the Lagrangian with the covariant derivative in terms of the SU (2) L and U (1) Y gauge fields W a µ and B µ , respectively, with their corresponding couplings g and g . The bilinear terms in the fermion fields lead to mass matrices for the 2/3, −1/3 and 5/3 charged fermions, when the Higgs field is shifted by its VEV H , H = H + h. The mass matrices can be diagonalized by means of a bi-unitary transformation. The 2-fermion couplings to one and two Higgs bosons are obtained by expanding the mass matrices in the interaction eigenstates up to first, respectively, second order in the Higgs field, and subsequent transformation into the mass eigenstate basis. The mass matrices and the coupling matrices of one Higgs boson to two bottom-like and top-like states can be found in Ref. [102]. In the Appendix A we give the coupling matrices for the 2-Higgs-2-fermion couplings and, for completeness, repeat the matrices given in Ref. [102].

Next-to-leading Order QCD Corrections to Higgs Pair Production in Composite Higgs Models
The NLO QCD corrections to Higgs pair production in the SM have been computed in Ref. [4] by applying the heavy top approximation, in which the heavy fermion loops are replaced by effective vertices of gluons to Higgs bosons. These can be obtained by means of the low-energy theorem (LET) [103,104]. The Higgs field is treated here as a background field, and the field-dependent mass of each heavy particle is taken into account in the gluon self-interactions at higher orders.
The LET provides the zeroth order in an expansion in small external momenta. Since in Higgs pair production the requirements for such an expansion are not fulfilled sufficiently reliably, it fails to give accurate results for the cross section at LO [81]. In the context of composite Higgs models, the discrepancy between the LO cross section with full top quark mass dependence and the LO cross section in the LET approximation is even worse [74]. For relative higher order corrections the LET approximation should, however, become better, if the LO order cross section is taken into account with full mass dependence. This is because the dominant corrections given by the soft and collinear gluon corrections factorize from the LO cross section generating a part independent of the masses of the heavy loop particles relative to the LO cross section. This was confirmed in Ref. [76] by including higher terms in the expansion of the cross section in small external momenta. Based on these findings, in this section we will give the NLO QCD corrections for Higgs pair production in composite Higgs models in the LET approach.
The expression of the LO gluon fusion into Higgs pairs in a composite Higgs model with heavy top partners has been given in [74]. It can be taken over here, by simply extending the sum to include also the bottom quark and its partners. We summarize here the most important features and refer to [74] for more details. The generic diagrams that contribute to the process at LO are depicted in Fig. 1. Besides the new 2-Higgs-2-fermion coupling f f hh the additional top and bottom partners in the loops have to be taken into account. These lead also to new box diagrams involving off-diagonal Yukawa couplings, with, respectively, the top and its heavy charge-2/3 partners or the bottom and its heavy partners of charge −1/3. The hadronic cross section is obtained by convolution with the parton distribution functions f g of the gluon in the proton, where s denotes the squared hadronic c.m. energy, µ F the factorization scale and in terms of the Higgs boson mass m h . The partonic LO cross section can be cast into the form with the strong coupling constant α s at the renormalization scale µ R . We have introduced the Mandelstam variableŝ in terms of the scattering angle θ in the partonic c.m. system with the invariant Higgs pair mass Q and the relative velocity The integration limits at cos θ = ±1 are given bŷ The form factors read The triangle and box form factors F ∆ , F , F ,5 , G and G ,5 can be found in the appendices of [74,105]. 7 The sum runs up to n t = 5 for the top quark and its charge-2/3 partners and up to n b = 4 in the bottom sector. The couplings are defined as and 7 The form factors F∆, F and G relate to those given in Ref. [82] for the SM case as where G hqq,ij and G hhqq,ij denote the (ith,jth) matrix elements of the coupling matrices in Eq. (58) of the appendix. The triangle factor C i,∆ reads in the MCHM10 as given in the MCHM5. In the SM and in the composite Higgs models MCHM4 and MCHM5 involving solely the Higgs non-linearities and no heavy fermionic resonances, no sum over heavy top and bottom partners contributes and only a sum over the top and bottom running in the loop has to be performed, i.e. n t = n b = 1, with m i = m j = m q and q = t, b, and hence also g hq i q j ,5 = 0 for SM, MCHM4 and MCHM5.
The Yukawa couplings read and for the 2-Higgs-2-fermion coupling we have while the Higgs self-coupling becomes The Feynman diagrams contributing to Higgs pair production at NLO QCD are shown in Fig. 2. The blob in the figure marks the effective vertices of gluons to Higgs boson(s). The first three Feynman diagrams show the virtual contributions. The remaining Feynman diagrams of Fig. 2 display the real corrections generically ordered by the initial states gg, gq and qq. At NLO the cross section is then given by σ NLO (pp → hh + X) = σ LO + ∆σ virt + ∆σ gg + ∆σ gq + ∆σ qq .
The individual contributions in Eq. (30) read with the Altarelli-Parisi splitting functions given by [106] P gg (z) = 6 and N F = 5 in our case. The real corrections ∆σ gg , ∆σ gq and ∆σ qq have straightforwardly been obtained from Ref. [4] by replacing the LO cross section of the SM with the LO cross section for composite Higgs models. The calculation of ∆σ virt is a bit more involved. While the first two diagrams factorize from the LO cross section and can hence directly be taken over from the SM, the third diagram in Fig. 2 does not factorize and needs to be recalculated for the composite Higgs case. The virtual coefficient C is then found to be with The first line in Eq. (37) corresponds to the NLO contribution from the first two diagrams in Fig. 2, while the second line corresponds to the NLO contribution from the third diagram of Fig. 2. The factor (g eff hgg ) 2 stems from the two effective Higgs-gluon-gluon vertices in diagram 3 of Fig. 2. This vertex is obtained by integrating out all heavy loop particles in the loop-induced Higgs coupling to gluons defined in Eq. (6) with g hgg ≡ g eff hgg and The first term is the sum over the normalized top quark and top partner couplings and the second term the sum over the normalized bottom partner couplings to the Higgs boson, excluding consistently the light bottom quark contribution from the loop. The composite Higgs cross sections for MCHM4, MCHM5 and for the composite Higgs model with heavy top and bottom partners, including the NLO corrections have been implemented in HPAIR [107]. In order to exemplify the impact of the NLO QCD corrections, we consider the simple case with the pure Higgs non-linearities only and the fermions transforming in the fundamental representation, i.e. the benchmark model MCHM5, see Table 1. The coupling g eff hgg then reduces to g MCHM5 hgg = (1 − 2ξ)/ √ 1 − ξ and the remaining couplings are given in Eqs. (26)- (29). We define the K-factors for the total cross section and the individual contributions as The cross section at LO is computed with the full quark mass dependences. As the NLO cross section in the LET approximation only includes top quark contributions 8 The renormalization and factorization scales are set to µ R = µ F = m hh /2, where m hh denotes the invariant Higgs pair mass.  contributions. As can be inferred from the plot, the real and virtual corrections of the gg initial state make up the bulk of the QCD corrections. The qg and the qq initiated real radiation diagrams only lead to a small correction. The K-factor is almost independent of ξ. In the real corrections, the Born cross section, which shows the only dependence on ξ, almost completely drops out numerically. For the virtual contributions some dependence on ξ may be expected. The virtual correction due to the constant term in C, i.e. the first line in Eq. (37) does not develop any dependence on ξ, however, as it factorizes from the LO cross section. The dependence of ξ can only emerge from the second line in Eq. (37), which, however, is numerically suppressed. This is already the case in the SM, where the corresponding term contributes less than 3% to ∆σ virt . This result also holds true for the case were the heavy quark partners are explicitly included. In composite Higgs models, the NLO QCD corrections to Higgs pair production can hence well be approximated by multiplying the full LO cross section of the composite Higgs model under consideration with the SM K-factor. Figure 3 can also be obtained by using the results of Ref. [89]. Note however, that the effects of heavy top and bottom partners in the effective field theory computation of Ref. [89] have to be added to the top quark contribution, encoded into the Wilson coefficients in front of the operators hG µν G µν and hhG µν G µν .

Numerical Analysis of New Physics Effects in Higgs Pair Production via Gluon Fusion
Having derived the NLO QCD corrections, we can now turn to the analysis of NP effects in Higgs pair production. We assume that no NP is found before Higgs pair production becomes accessible. This means that we require deviations in the Higgs boson couplings with respect to the SM to be smaller than the projected sensitivities of the coupling measurements at an integrated luminosity of 300 fb −1 and 3000 fb −1 , respectively. For the projected sensitivities we take the numbers reported in Ref. [109]. Similar numbers can be found in Refs. [110]. In our analysis we focus on the most promising final states, given by bbγγ and bbτ + τ − [8][9][10]13].
We call Higgs pair production to be sensitive to NP if the difference between the number of signal events S in the considered NP model and the corresponding number S SM in the SM exceeds a minimum of 3 statistical standard deviations, i.e.
with β = 3 for a 3σ deviation. The signal events are obtained as where BR denotes the branching ratio into the respective final states, L the integrated luminosity and A the acceptance due to cuts applied on the cross section. The acceptances have been extracted from Ref. [13] for the bbγγ and bbτ + τ − final states. The acceptance for the BSM signal only changes slightly, as we explicitly checked.
In specific models, the correlations of the couplings will lead to stronger bounds on the parameters. In particular in the MCHM4 and MCHM5 as introduced in section 2, the only new parameter is ξ. The value of ξ can hence strongly be restricted by Higgs coupling measurements [111]. Based on these estimates, we give in Table 2 the maximal values for the cross section times branching ratio. In the fourth and sixth columns we report whether the process within MCHM4, respectively, MCHM5 can be distinguished from the SM cross section by more than 3σ according to the criteria given in Eq. (42) for β = 3. In the check of Eq. (42) we took into account the slight change in the acceptance of the signal rate for the composite Higgs models. Due to the coupling modifications and the new diagram from the 2-Higgs-2-fermion coupling the applied cuts in the analysis of Ref. [13] affect the cross section in a slightly different way.
The table shows, that with the projected precision on ξ at high luminosities Higgs pair production in both MCHM4 and MCHM5 leads to cross sections too close to the SM value to be distinguishable from the SM case. Although with the present bounds on ξ Higgs pair production in MCHM5 differs by more than 3σ from the SM prediction, the corresponding cross section is too small to be measurable, so that first signs of NP through this process are precluded.

Numerical Analysis for MCHM10
We consider the MCHM10 with one multiplet of fermionic resonances below the cut-off. In this model, with more than one parameter determining the Higgs coupling modifications, there is more freedom and a larger allowed parameter space (see Ref. [102] for a thorough analysis). This implies,  Table 2: Values of the cross section times branching ratio in the MCHM4 and MCHM5 for the maximal allowed values of ξ at 95% C.L.
[112] and for the projected values at L = 300 fb −1 and L = 3000 fb −1 of Ref. [109]. The fourth and sixth columns decide whether the Higgs production cross section will develop a deviation to the SM Higgs pair production cross section of more than 3σ according to the criteria of Eq. (42).
that the sensitivity on the Higgs couplings is less constrained. The numbers of the projected sensitivities are taken from Table I in Ref. [109]. Additionally, we need to take into account the bounds from the direct searches for new fermions. Currently, exotic new fermions with charge 5/3 are excluded up to masses of m χ ≤ 840 GeV [113], bottom partners up to masses of m B ≤ 900 GeV [114] and top partners with masses of m T ≤ 950 GeV [115]. Note that the latter two limits on the masses depend on the branching ratios of the bottom and top partner, respectively. These limits are based on pair production of the new heavy fermions. First studies for single production of a new vector-like fermion were performed in Refs. [116] and can potentially be more important at large energies [117] but are more model-dependent. Due to this model dependence it is difficult to estimate the LHC reach on single production for our case. Hence we will only use the estimated reach on new vector-like fermions in pair production. In Refs. [118,119] the potential reach of the LHC for charged-2/3 fermions, depending on their branching ratios is estimated. Following [119] we use the reach m T 1.3 TeV for L = 300 fb −1 and m T 1.5 TeV for L = 3000 fb −1 . The potential reach for bottom partners is m B 1 TeV for L = 300 fb −1 and m B 1.5 TeV for L = 3000 fb −1 [120]. We estimate the additional sensitivity for the reach of exotic new fermions by multiplying the excluded cross section at √ s = 8 TeV with [121] r = σ BKG (14 TeV) σ BKG (8 TeV) where L LHC8 and L LHC14 denote the integrated luminosities of the LHC run at √ s = 8 and 14 TeV, respectively. This implies a reach of m χ ≈ 1370 GeV at L = 300 fb −1 and m χ ≈ 1550 GeV at L = 3000 fb −1 . For the background estimate we only considered the dominant background ttW ± [122]. The background cross section was computed with MadGraph5 [123]. Although the assumption of stronger projections on the reach of new fermion masses of up to 2 TeV [124] will lead to a reduced number of points allowed by the constraints we are imposing, it will not change our final conclusion, as we checked explicitly. Note also that in composite Higgs models there is a connection between the Higgs boson mass and the fermionic resonances [125,126]. Reference [126] finds that the mass of the lightest top partner m T lightest should be lighter than with N c = 3 denoting the number of colors. This bound automatically eliminates large values of ξ. In our analysis we allow for finetuning, hence small values of ξ, and we will not apply this bound. For the analysis we performed a scan over the parameter space of the model by varying the parameters in the range 9 We excluded points that do not fulfill |V tb | > 0.92 [127] and the electroweak precision tests (EWPTs) at 99% C.L. using the results of Ref. [102].
In Fig. 4 we show the NLO Higgs pair production cross section via gluon fusion as a function of ξ. The color code in the plots indicates whether the points are distinguishable from the SM according to the criteria given in Eq. (42), with the blue points being distinguishable and the grey points not. The upper plots are for the bbτ + τ − final state, the lower plots for the bbγγ final state, for L = 300 fb −1 (left) and L = 3000 fb −1 (right), respectively. The upper branch in the plots corresponds to the parameters y < 0 and 0 < R < 1 with R = (M 10 + f y/2)/M 10 . This means that at LO of the mass matrix expansion in v/f , the lightest fermion partner originates from the SU (2) bi-doublet. The lower branch corresponds to the cases y < 0 and R < 0 as well as y > 0 implying R > 1.
The plots only show points for which we cannot see new physics anywhere else meaning we require that their deviations in the Higgs couplings that can be tested at the LHC are smaller than the expected sensitivities and that the masses of the new fermionic resonances are above the estimated reach of direct searches. 10 The requirement for only small deviations in the Higgs couplings directly restricts the possible values of ξ to be smaller than 0.071 and 0.059 for L = 300 fb −1 and L = 3000 fb −1 , respectively. The value of ξ is restricted more strongly than in the MCHM4 due to the different coupling modifications, which, considering pure non-linearities, are for the Higgs-fermion couplings (1 − 2ξ)/ √ 1 − ξ in MCHM10 and √ 1 − ξ in the MCHM4. Although the interplay of the various additional parameters in MCHM10 allows for some tuning in the Higgs-bottom (and also Higgs-top) coupling, this is not the case for the Higgs-tau coupling. Comparing the MCHM10 with the MCHM5, the Higgs-fermion couplings are modified in the same way, barring the effects from the additional fermions. The increased number of parameters due to the heavy fermions, however, allows for more freedom to accommodate the data, so that here the constraint is weaker in the MCHM10. The plots show that at L = 300 fb −1 we cannot expect to discover NP for the first time in Higgs pair production in the bbγγ final state, while in the bbτ + τ − final state NP could show up for the first time in Higgs pair production. For L = 3000 fb −1 , we could find both in the bbτ + τ − and the bbγγ final state points which lead to large enough deviations from the SM case to be sensitive to NP for the first time in Higgs pair production. These results can be explained with the increased signal rate in the cases that are sensitive, as can be inferred from Fig. 5. The plots show for the parameter points displayed in Fig. 4 the corresponding number of signal events for Higgs pair production in the bbτ + τ − and bbγγ final states, respectively, after applying the acceptance of the cuts and multiplication with the two options of integrated luminosity. The blue points clearly deviate by more than 3σ from the SM curve.

Invariant Mass Distributions
Finally, in this section we discuss NP effects in invariant Higgs mass distributions. The measurement of distributions can give information on anomalous couplings [69] or the underlying ultraviolet source of NP [128]. Even though they are difficult to be measured due to the small numbers of signal events, they are important observables for the NP search. In the following we will show the impact of composite Higgs models on the distributions.  Note, however, that the shape of the invariant mass distributions hardly changes from LO to NLO, since in the LET approximation the LO cross section mainly factorizes from the NLO contributions, as discussed in section 3. The parameters have been chosen such that in the left plot we allow for a larger value of ξ, while the mass of the lightest top partner of m T,lightest = 5441 GeV is much larger than compared to the case shown in the right plot with m T,lightest = 1636 GeV. As can be inferred from these plots, the largest effect on the distributions originates from the pure non-linearities, i.e. the value of ξ, while the influence of the fermionic resonances on the shape of the invariant mass distributions is small. Note also that the main effect on the distributions emerges from the new tthh coupling.

Conclusions
We presented the NLO QCD corrections to Higgs pair production via gluon fusion in the large quark mass approximation in composite Higgs models with and without new heavy fermionic resonances below the UV cut-off. We found that the K-factor of ∼ 1.7 is basically independent of the value of ξ and of the details of the fermion spectrum, as the LO cross section dominantly factorizes. The K-factor can hence directly be taken over from the SM to a good approximation. The size of the absolute value of the cross section, however, sensitively depends on the Higgs non-linearities and on the fermion spectrum.
With the results of our NLO calculation, we furthermore addressed the question of whether NP could emerge for the first time in Higgs pair production, taking into account the constraints on the Higgs couplings to SM particles and from direct searches for new heavy fermions. We focused on composite Higgs models and found that in simple models where only the Higgs non-linearities are considered, we cannot expect to be sensitive to NP for the first time in Higgs pair production. In models with a multiplet of new fermions below the cut-off, it turned out that there are regions where NP could indeed be seen for the first time in Higgs pair production. The subsequent study of the NLO invariant mass distributions demonstrated, that while there is some sensitivity to the Higgs non-linearities mainly due to the new 2-Higgs-2-fermion coupling, the effect of the heavy fermions on the shape of the distributions is much weaker. By applying optimized cuts the sensitivity to new physics effects may possibly be increased in future. the terms bilinear in the quark fields of Eq. (12) read and where U (t/b/χ) L,R denote the transformations that diagonalize the mass matrix in the top, bottom and charge-5/3 (χ) sector, respectively. Expansion of the mass matrices Eqs. (50)- (52) in the interaction eigenstates up to first order in the Higgs field leads to the Higgs coupling matricesG htt to a pair of charge-2/3 quarks andG hbb to a quark pair of charge −1/3, respectively, given by Expansion up to second order in the Higgs field yields the 2-Higgs-2-fermion coupling matrices G hhtt andG hhbb , that can be cast into the form The coupling matrices in the mass eigenstate basis are obtained by rotation with the unitary matrices defined in Eq. (53), i.e. (q = t, b)